
  module cldwat2m_macro

  !--------------------------------------------------- !
  ! Purpose     : CAM Interface for Cloud Macrophysics !
  ! Author      : Sungsu Park                          !
  ! Description : Park et al. 2010.                    !
  ! For questions, contact Sungsu Park                 !
  !                        e-mail : sungsup@ucar.edu   !
  !                        phone  : 303-497-1375       !
  !--------------------------------------------------- !

   use shr_kind_mod,     only: r8=>shr_kind_r8
   use spmd_utils,       only: masterproc
   use ppgrid,           only: pcols, pver, pverp
   use abortutils,       only: endrun
   use wv_saturation,    only: estblf, hlatv, tmin, hlatf, rgasv, pcf, cp, epsqs, ttrice, & 
                               vqsatd2_water, vqsatd2_water_single, polysvp
   use cam_history,      only: addfld, add_default, phys_decomp, outfld 
   use cam_logfile,      only: iulog

   implicit none
   private
   save

   public :: ini_macro, mmacro_pcond

   ! -------------- !
   ! Set Parameters !
   ! -------------- !

   ! -------------------------- !
   ! Parameters for Ice Stratus !
   ! -------------------------- !

   integer,  parameter          :: iceopt       = 5             ! Option for ice cloud closure. 5 : Modified Slingo Formula. For other options, modify each 2 subroutine.
   real(r8), public,  parameter :: rhmini       = 0.80_r8       ! Minimum rh for ice cloud fraction > 0.
   real(r8), public,  parameter :: rhmaxi       = 1.1_r8        ! rhi at which ice cloud fraction = 1.
   real(r8), parameter          :: qist_min     = 1.e-7_r8      ! Minimum in-stratus ice IWC constraint [ kg/kg ]
   real(r8), parameter          :: qist_max     = 5.e-3_r8      ! Maximum in-stratus ice IWC constraint [ kg/kg ]

   ! ----------------------------- !
   ! Parameters for Liquid Stratus !
   ! ----------------------------- !

   logical,  parameter          :: CAMstfrac    = .false.       ! If .true. (.false.), use Slingo (triangular PDF-based) liquid stratus fraction
   logical,  parameter          :: freeze_dry   = .false.       ! If .true., use 'freeze dry' in liquid stratus fraction formula
   real(r8), parameter          :: qlst_min     = 2.e-5_r8      ! Minimum in-stratus LWC constraint [ kg/kg ]
   real(r8), parameter          :: qlst_max     = 3.e-3_r8      ! Maximum in-stratus LWC constraint [ kg/kg ]
   real(r8), parameter          :: cc           = 0.1_r8        ! For newly formed/dissipated in-stratus CWC ( 0 <= cc <= 1 )
   real(r8), parameter          :: premib       = 700.e2_r8     ! Bottom height for mid-level liquid stratus fraction
   integer,  parameter          :: niter        = 2             ! For iterative computation of QQ with 'ramda' below.
   real(r8), parameter          :: ramda        = 0.5_r8        ! Explicit : ramda = 0, Implicit : ramda = 1 ( 0<= ramda <= 1 )
   real(r8), private            :: rhminl                       ! Critical RH for low-level  liquid stratus clouds
   real(r8), private            :: rhminh                       ! Critical RH for high-level liquid stratus clouds
   real(r8), private            :: premit                       ! Top    height for mid-level liquid stratus fraction

   contains

   ! -------------- !
   ! Initialization !
   ! -------------- !

   subroutine ini_macro

   !--------------------------------------------------------------------- !
   !                                                                      ! 
   ! Purpose: Initialize constants for the liquid stratiform macrophysics !
   !          called from stratiform.F90.                                 !  
   ! Author:  Sungsu Park, Dec.01.2009.                                   !
   !                                                                      !
   !--------------------------------------------------------------------- !

   use cloud_fraction, only: cldfrc_getparams

   call cldfrc_getparams(rhminl_out = rhminl , rhminh_out = rhminh, premit_out = premit)

   if( masterproc ) then
       write(iulog,*) 'ini_macro: tuning parameters : rhminl = ', rhminl, &
                                                     'rhminh = ', rhminh, & 
                                                     'premit = ', premit, & 
                                                     'premib = ', premib
   end if

   return
   end subroutine ini_macro

   ! ------------------------------ !
   ! Stratiform Liquid Macrophysics !
   ! ------------------------------ !

   ! In the version, 'macro --> micro --> advective forcing --> macro...'
   ! A_...: only 'advective forcing' without 'microphysical tendency'
   ! C_...: only 'microphysical tendency'
   ! D_...: only 'detrainment of cumulus condensate'  
   ! So, 'A' and 'C' are exclusive. 

   subroutine mmacro_pcond( lchnk      , ncol       , dt         , p            , dp         ,              &
                            T0         , qv0        , ql0        , qi0          , nl0        , ni0        , &
                            A_T        , A_qv       , A_ql       , A_qi         , A_nl       , A_ni       , &
                            C_T        , C_qv       , C_ql       , C_qi         , C_nl       , C_ni       , C_qlst     , &
                            D_T        , D_qv       , D_ql       , D_qi         , D_nl       , D_ni       , &
                            a_cud      , a_cu0      , landfrac   , snowh        ,                           & 
                            s_tendout  , qv_tendout , ql_tendout , qi_tendout   , nl_tendout , ni_tendout , &
                            qme        , qvadj      , qladj      , qiadj        , qllim      , qilim      , &
                            cld        , al_st_star , ai_st_star , ql_st_star   , qi_st_star )

   use constituents,     only : qmin
   use time_manager,     only : is_first_step, get_nstep
   use phys_debug_util,  only : phys_debug_col

   implicit none

   integer   icol
   integer,  intent(in)    :: lchnk                        ! Chunk number
   integer,  intent(in)    :: ncol                         ! Number of active columns

   ! Input-Output variables

   real(r8), intent(inout) :: T0(pcols,pver)               ! Temperature [K]
   real(r8), intent(inout) :: qv0(pcols,pver)              ! Grid-mean water vapor specific humidity [kg/kg]
   real(r8), intent(inout) :: ql0(pcols,pver)              ! Grid-mean liquid water content [kg/kg]
   real(r8), intent(inout) :: qi0(pcols,pver)              ! Grid-mean ice water content [kg/kg]
   real(r8), intent(inout) :: nl0(pcols,pver)              ! Grid-mean number concentration of cloud liquid droplet [#/kg]
   real(r8), intent(inout) :: ni0(pcols,pver)              ! Grid-mean number concentration of cloud ice    droplet [#/kg]

   ! Input variables

   real(r8), intent(in)    :: dt                           ! Model integration time step [s]
   real(r8), intent(in)    :: p(pcols,pver)                ! Pressure at the layer mid-point [Pa]
   real(r8), intent(in)    :: dp(pcols,pver)               ! Pressure thickness [Pa] > 0

   real(r8), intent(in)    :: A_T(pcols,pver)              ! Non-microphysical advective external forcing of T  [K/s]
   real(r8), intent(in)    :: A_qv(pcols,pver)             ! Non-microphysical advective external forcing of qv [kg/kg/s]
   real(r8), intent(in)    :: A_ql(pcols,pver)             ! Non-microphysical advective external forcing of ql [kg/kg/s]
   real(r8), intent(in)    :: A_qi(pcols,pver)             ! Non-microphysical advective external forcing of qi [kg/kg/s]
   real(r8), intent(in)    :: A_nl(pcols,pver)             ! Non-microphysical advective external forcing of nl [#/kg/s]
   real(r8), intent(in)    :: A_ni(pcols,pver)             ! Non-microphysical advective external forcing of ni [#/kg/s] 

   real(r8), intent(in)    :: C_T(pcols,pver)              ! Microphysical advective external forcing of T  [K/s]
   real(r8), intent(in)    :: C_qv(pcols,pver)             ! Microphysical advective external forcing of qv [kg/kg/s]
   real(r8), intent(in)    :: C_ql(pcols,pver)             ! Microphysical advective external forcing of ql [kg/kg/s]
   real(r8), intent(in)    :: C_qi(pcols,pver)             ! Microphysical advective external forcing of qi [kg/kg/s]
   real(r8), intent(in)    :: C_nl(pcols,pver)             ! Microphysical advective external forcing of nl [#/kg/s]
   real(r8), intent(in)    :: C_ni(pcols,pver)             ! Microphysical advective external forcing of ni [#/kg/s] 
   real(r8), intent(in)    :: C_qlst(pcols,pver)           ! Microphysical advective external forcing of ql within liquid stratus [kg/kg/s]

   real(r8), intent(in)    :: D_T(pcols,pver)              ! Cumulus detrainment external forcing of T  [K/s]
   real(r8), intent(in)    :: D_qv(pcols,pver)             ! Cumulus detrainment external forcing of qv [kg/kg/s]
   real(r8), intent(in)    :: D_ql(pcols,pver)             ! Cumulus detrainment external forcing of ql [kg/kg/s]
   real(r8), intent(in)    :: D_qi(pcols,pver)             ! Cumulus detrainment external forcing of qi [kg/kg/s]
   real(r8), intent(in)    :: D_nl(pcols,pver)             ! Cumulus detrainment external forcing of nl [#/kg/s]
   real(r8), intent(in)    :: D_ni(pcols,pver)             ! Cumulus detrainment external forcing of qi [#/kg/s] 

   real(r8), intent(in)    :: a_cud(pcols,pver)            ! Old cumulus fraction before update
   real(r8), intent(in)    :: a_cu0(pcols,pver)            ! New cumulus fraction after update

   real(r8), intent(in)    :: landfrac(pcols)              ! Land fraction
   real(r8), intent(in)    :: snowh(pcols)                 ! Snow depth (liquid water equivalent)

   ! Output variables

   real(r8), intent(out)   :: s_tendout(pcols,pver)        ! Net tendency of grid-mean s  from 'Micro+Macro' processes [J/kg/s]
   real(r8), intent(out)   :: qv_tendout(pcols,pver)       ! Net tendency of grid-mean qv from 'Micro+Macro' processes [kg/kg/s]
   real(r8), intent(out)   :: ql_tendout(pcols,pver)       ! Net tendency of grid-mean ql from 'Micro+Macro' processes [kg/kg/s]
   real(r8), intent(out)   :: qi_tendout(pcols,pver)       ! Net tendency of grid-mean qi from 'Micro+Macro' processes [kg/kg/s]
   real(r8), intent(out)   :: nl_tendout(pcols,pver)       ! Net tendency of grid-mean nl from 'Micro+Macro' processes [#/kg/s]
   real(r8), intent(out)   :: ni_tendout(pcols,pver)       ! Net tendency of grid-mean ni from 'Micro+Macro' processes [#/kg/s]

   real(r8), intent(out)   :: qme  (pcols,pver)            ! Net condensation rate [kg/kg/s]
   real(r8), intent(out)   :: qvadj(pcols,pver)            ! adjustment tendency from "positive_moisture" call (vapor)
   real(r8), intent(out)   :: qladj(pcols,pver)            ! adjustment tendency from "positive_moisture" call (liquid)
   real(r8), intent(out)   :: qiadj(pcols,pver)            ! adjustment tendency from "positive_moisture" call (ice)
   real(r8), intent(out)   :: qllim(pcols,pver)            ! tendency from "instratus_condensate" call (liquid)
   real(r8), intent(out)   :: qilim(pcols,pver)            ! tendency from "instratus_condensate" call (ice)

   real(r8), intent(out)   :: cld(pcols,pver)              ! Net cloud fraction ( 0 <= cld <= 1 )
   real(r8), intent(out)   :: al_st_star(pcols,pver)       ! Physical liquid stratus fraction
   real(r8), intent(out)   :: ai_st_star(pcols,pver)       ! Physical ice stratus fraction
   real(r8), intent(out)   :: ql_st_star(pcols,pver)       ! In-stratus LWC [kg/kg] 
   real(r8), intent(out)   :: qi_st_star(pcols,pver)       ! In-stratus IWC [kg/kg] 

   ! --------------- !
   ! Local variables !
   ! --------------- !
 
   integer :: i, j, k, iter, ii, jj                        ! Loop indexes

   ! Thermodynamic state variables

   real(r8) T(pcols,pver)                                  ! Temperature of equilibrium reference state from which 'Micro & Macro' are computed [K]
   real(r8) T1(pcols,pver)                                 ! Temperature after 'fice_force' on T01  
   real(r8) T_0(pcols,pver)                                ! Temperature after 'instratus_condensate' on T1
   real(r8) T_05(pcols,pver)                               ! Temperature after 'advection' on T_0 
   real(r8) T_prime0(pcols,pver)                           ! Temperature after 'Macrophysics (QQ)' on T_05star
   real(r8) T_dprime(pcols,pver)                           ! Temperature after 'fice_force' on T_prime
   real(r8) T_star(pcols,pver)                             ! Temperature after 'instratus_condensate' on T_dprime

   real(r8) qv(pcols,pver)                                 ! Grid-mean qv of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
   real(r8) qv1(pcols,pver)                                ! Grid-mean qv after 'fice_force' on qv01  
   real(r8) qv_0(pcols,pver)                               ! Grid-mean qv after 'instratus_condensate' on qv1
   real(r8) qv_05(pcols,pver)                              ! Grid-mean qv after 'advection' on qv_0 
   real(r8) qv_prime0(pcols,pver)                          ! Grid-mean qv after 'Macrophysics (QQ)' on qv_05star
   real(r8) qv_dprime(pcols,pver)                          ! Grid-mean qv after 'fice_force' on qv_prime
   real(r8) qv_star(pcols,pver)                            ! Grid-mean qv after 'instratus_condensate' on qv_dprime

   real(r8) ql(pcols,pver)                                 ! Grid-mean ql of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
   real(r8) ql1(pcols,pver)                                ! Grid-mean ql after 'fice_force' on ql01  
   real(r8) ql_0(pcols,pver)                               ! Grid-mean ql after 'instratus_condensate' on ql1
   real(r8) ql_05(pcols,pver)                              ! Grid-mean ql after 'advection' on ql_0 
   real(r8) ql_prime0(pcols,pver)                          ! Grid-mean ql after 'Macrophysics (QQ)' on ql_05star
   real(r8) ql_dprime(pcols,pver)                          ! Grid-mean ql after 'fice_force' on ql_prime
   real(r8) ql_star(pcols,pver)                            ! Grid-mean ql after 'instratus_condensate' on ql_dprime

   real(r8) qi(pcols,pver)                                 ! Grid-mean qi of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
   real(r8) qi1(pcols,pver)                                ! Grid-mean qi after 'fice_force' on qi01  
   real(r8) qi_0(pcols,pver)                               ! Grid-mean qi after 'instratus_condensate' on qi1
   real(r8) qi_05(pcols,pver)                              ! Grid-mean qi after 'advection' on qi_0 
   real(r8) qi_prime0(pcols,pver)                          ! Grid-mean qi after 'Macrophysics (QQ)' on qi_05star
   real(r8) qi_dprime(pcols,pver)                          ! Grid-mean qi after 'fice_force' on qi_prime
   real(r8) qi_star(pcols,pver)                            ! Grid-mean qi after 'instratus_condensate' on qi_dprime

   real(r8) nl(pcols,pver)                                 ! Grid-mean nl of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
   real(r8) nl1(pcols,pver)                                ! Grid-mean nl after 'fice_force' on nl01  
   real(r8) nl_0(pcols,pver)                               ! Grid-mean nl after 'instratus_condensate' on nl1
   real(r8) nl_05(pcols,pver)                              ! Grid-mean nl after 'advection' on nl_0 
   real(r8) nl_prime0(pcols,pver)                          ! Grid-mean nl after 'Macrophysics (QQ)' on nl_05star
   real(r8) nl_dprime(pcols,pver)                          ! Grid-mean nl after 'fice_force' on nl_prime
   real(r8) nl_star(pcols,pver)                            ! Grid-mean nl after 'instratus_condensate' on nl_dprime

   real(r8) ni(pcols,pver)                                 ! Grid-mean ni of equilibrium reference state from which 'Micro & Macro' are computed [kg/kg]
   real(r8) ni1(pcols,pver)                                ! Grid-mean ni after 'fice_force' on ni01  
   real(r8) ni_0(pcols,pver)                               ! Grid-mean ni after 'instratus_condensate' on ni1
   real(r8) ni_05(pcols,pver)                              ! Grid-mean ni after 'advection' on ni_0 
   real(r8) ni_prime0(pcols,pver)                          ! Grid-mean ni after 'Macrophysics (QQ)' on ni_05star
   real(r8) ni_dprime(pcols,pver)                          ! Grid-mean ni after 'fice_force' on ni_prime
   real(r8) ni_star(pcols,pver)                            ! Grid-mean ni after 'instratus_condensate' on ni_dprime

   real(r8) a_st(pcols,pver)                               ! Stratus fraction of equilibrium reference state 
   real(r8) a_st_0(pcols,pver)                             ! Stratus fraction at '_0' state
   real(r8) a_st_star(pcols,pver)                          ! Stratus fraction at '_star' state

   real(r8) al_st(pcols,pver)                              ! Liquid stratus fraction of equilibrium reference state 
   real(r8) al_st_0(pcols,pver)                            ! Liquid stratus fraction at '_0' state
   real(r8) al_st_nc(pcols,pver)                           ! Non-physical liquid stratus fraction in the non-cumulus pixels

   real(r8) ai_st(pcols,pver)                              ! Ice stratus fraction of equilibrium reference state 
   real(r8) ai_st_0(pcols,pver)                            ! Ice stratus fraction at '_0' state
   real(r8) ai_st_nc(pcols,pver)                           ! Non-physical ice stratus fraction in the non-cumulus pixels

   real(r8) ql_st(pcols,pver)                              ! In-stratus LWC of equilibrium reference state [kg/kg] 
   real(r8) ql_st_0(pcols,pver)                            ! In-stratus LWC at '_0' state

   real(r8) qi_st(pcols,pver)                              ! In-stratus IWC of equilibrium reference state [kg/kg] 
   real(r8) qi_st_0(pcols,pver)                            ! In-stratus IWC at '_0' state

 ! Cumulus properties 

   real(r8) dacudt(pcols,pver)
   real(r8) a_cu(pcols,pver)

 ! Adjustment tendency in association with 'positive_moisture'

   real(r8) Tten_pwi1(pcols,pver)                          ! Pre-process T  tendency of input equilibrium state [K/s] 
   real(r8) qvten_pwi1(pcols,pver)                         ! Pre-process qv tendency of input equilibrium state [kg/kg/s]
   real(r8) qlten_pwi1(pcols,pver)                         ! Pre-process ql tendency of input equilibrium state [kg/kg/s]
   real(r8) qiten_pwi1(pcols,pver)                         ! Pre-process qi tendency of input equilibrium state [kg/kg/s]
   real(r8) nlten_pwi1(pcols,pver)                         ! Pre-process nl tendency of input equilibrium state [#/kg/s]
   real(r8) niten_pwi1(pcols,pver)                         ! Pre-process ni tendency of input equilibrium state [#/kg/s] 

   real(r8) Tten_pwi2(pcols,pver)                          ! Post-process T  tendency of provisional equilibrium state [K/s] 
   real(r8) qvten_pwi2(pcols,pver)                         ! Post-process qv tendency of provisional equilibrium state [kg/kg/s]
   real(r8) qlten_pwi2(pcols,pver)                         ! Post-process ql tendency of provisional equilibrium state [kg/kg/s]
   real(r8) qiten_pwi2(pcols,pver)                         ! Post-process qi tendency of provisional equilibrium state [kg/kg/s]
   real(r8) nlten_pwi2(pcols,pver)                         ! Post-process nl tendency of provisoonal equilibrium state [#/kg/s]
   real(r8) niten_pwi2(pcols,pver)                         ! Post-process ni tendency of provisional equilibrium state [#/kg/s] 

   real(r8) A_T_adj(pcols,pver)                            ! After applying external advective forcing [K/s]
   real(r8) A_qv_adj(pcols,pver)                           ! After applying external advective forcing [kg/kg/s]
   real(r8) A_ql_adj(pcols,pver)                           ! After applying external advective forcing [kg/kg/s]
   real(r8) A_qi_adj(pcols,pver)                           ! After applying external advective forcing [kg/kg/s]
   real(r8) A_nl_adj(pcols,pver)                           ! After applying external advective forcing [#/kg/s]
   real(r8) A_ni_adj(pcols,pver)                           ! After applying external advective forcing [#/kg/s]

 ! Adjustment tendency in association with 'instratus_condensate'

   real(r8) QQw1(pcols,pver)                               ! Effective adjustive condensation into water due to 'instratus_condensate' [kg/kg/s]
   real(r8) QQi1(pcols,pver)                               ! Effective adjustive condensation into ice   due to 'instratus_condensate' [kg/kg/s]
   real(r8) QQw2(pcols,pver)                               ! Effective adjustive condensation into water due to 'instratus_condensate' [kg/kg/s]
   real(r8) QQi2(pcols,pver)                               ! Effective adjustive condensation into ice   due to 'instratus_condensate' [kg/kg/s]

   real(r8) QQnl1(pcols,pver)                              ! Tendency of nl associated with QQw1 only when QQw1<0 (net evaporation) [#/kg/s]
   real(r8) QQni1(pcols,pver)                              ! Tendency of ni associated with QQi1 only when QQw1<0 (net evaporation) [#/kg/s]
   real(r8) QQnl2(pcols,pver)                              ! Tendency of nl associated with QQw2 only when QQw2<0 (net evaporation) [#/kg/s]
   real(r8) QQni2(pcols,pver)                              ! Tendency of ni associated with QQi2 only when QQw2<0 (net evaporation) [#/kg/s]

 ! Macrophysical process tendency variables

   real(r8) QQ(pcols,pver)                                 ! Net condensation rate into water+ice           [kg/kg/s] 
   real(r8) QQw(pcols,pver)                                ! Net condensation rate into water               [kg/kg/s] 
   real(r8) QQi(pcols,pver)                                ! Net condensation rate into ice                 [kg/kg/s]
   real(r8) QQnl(pcols,pver)                               ! Tendency of nl associated with QQw both for condensation and evaporation [#/kg/s]
   real(r8) QQni(pcols,pver)                               ! Tendency of ni associated with QQi both for condensation and evaporation [#/kg/s]
   real(r8) ACnl(pcols,pver)                               ! Cloud liquid droplet (nl) activation tendency [#/kg/s]
   real(r8) ACni(pcols,pver)                               ! Cloud ice    droplet (ni) activation tendency [#/kg/s]

   real(r8) QQw_prev(pcols,pver)   
   real(r8) QQi_prev(pcols,pver)   
   real(r8) QQnl_prev(pcols,pver)  
   real(r8) QQni_prev(pcols,pver)  

   real(r8) QQw_prog(pcols,pver)   
   real(r8) QQi_prog(pcols,pver)   
   real(r8) QQnl_prog(pcols,pver)  
   real(r8) QQni_prog(pcols,pver)  

   real(r8) QQ_final(pcols,pver)                           
   real(r8) QQw_final(pcols,pver)                           
   real(r8) QQi_final(pcols,pver)                           
   real(r8) QQn_final(pcols,pver)                           
   real(r8) QQnl_final(pcols,pver)                          
   real(r8) QQni_final(pcols,pver)                          

   real(r8) QQ_all(pcols,pver)                             ! QQw_all    + QQi_all
   real(r8) QQw_all(pcols,pver)                            ! QQw_final  + QQw1  + QQw2  + qlten_pwi1 + qlten_pwi2 + A_ql_adj [kg/kg/s]
   real(r8) QQi_all(pcols,pver)                            ! QQi_final  + QQi1  + QQi2  + qiten_pwi1 + qiten_pwi2 + A_qi_adj [kg/kg/s]
   real(r8) QQn_all(pcols,pver)                            ! QQnl_all   + QQni_all
   real(r8) QQnl_all(pcols,pver)                           ! QQnl_final + QQnl1 + QQnl2 + nlten_pwi1 + nlten_pwi2 + ACnl [#/kg/s]
   real(r8) QQni_all(pcols,pver)                           ! QQni_final + QQni1 + QQni2 + niten_pwi1 + niten_pwi2 + ACni [#/kg/s]

 ! Coefficient for computing QQ and related processes

   real(r8) U(pcols,pver)                                  ! Grid-mean RH
   real(r8) U_nc(pcols,pver)                               ! Mean RH of non-cumulus pixels
   real(r8) G_nc(pcols,pver)                               ! d(U_nc)/d(a_st_nc)
   real(r8) F_nc(pcols,pver)                               ! A function of second parameter for a_st_nc
   real(r8) alpha                                          ! = 1/qsat
   real(r8) beta                                           ! = (qv/qsat**2)*dqsdT
   real(r8) betast                                         ! = alpha*dqsdT
   real(r8) gammal                                         ! = alpha + (hlatv/cp)*beta
   real(r8) gammai                                         ! = alpha + ((hlatv+hlatf)/cp)*beta
   real(r8) gammaQ                                         ! = alpha + (hlatv/cp)*beta
   real(r8) deltal                                         ! = 1 + a_st*(hlatv/cp)*(betast/alpha)
   real(r8) deltai                                         ! = 1 + a_st*((hlatv+hlatf)/cp)*(betast/alpha)
   real(r8) A_Tc                                           ! Advective external forcing of Tc [K/s]
   real(r8) A_qt                                           ! Advective external forcing of qt [kg/kg/s]
   real(r8) C_Tc                                           ! Microphysical forcing of Tc [K/s]
   real(r8) C_qt                                           ! Microphysical forcing of qt [kg/kg/s]
   real(r8) dTcdt                                          ! d(Tc)/dt      [K/s]
   real(r8) dqtdt                                          ! d(qt)/dt      [kg/kg/s]
   real(r8) dqtstldt                                       ! d(qt_alst)/dt [kg/kg/s]
   real(r8) dqidt                                          ! d(qi)/dt      [kg/kg/s]

   real(r8) dqlstdt                                        ! d(ql_st)/dt [kg/kg/s]
   real(r8) dalstdt                                        ! d(al_st)/dt  [1/s]
   real(r8) dastdt                                         ! d(a_st)/dt  [1/s]

   real(r8) anic                                           ! Fractional area of non-cumulus and non-ice stratus fraction
   real(r8) GG                                             ! G_nc(i,k)/anic

   real(r8) aa(2,2)
   real(r8) bb(2,1)

   real(r8) zeros(pcols,pver)

   real(r8) qmin1(pcols,pver)
   real(r8) qmin2(pcols,pver)
   real(r8) qmin3(pcols,pver)

   real(r8) esat_a(pcols)                                  ! Saturation water vapor pressure [Pa]
   real(r8) qsat_a(pcols)                                  ! Saturation water vapor specific humidity [kg/kg]
   real(r8) dqsdT_a(pcols)                                 ! dqsat/dT [kg/kg/K]
   real(r8) Twb_aw(pcols,pver)                             ! Wet-bulb temperature [K]
   real(r8) qvwb_aw(pcols,pver)                            ! Wet-bulb water vapor specific humidity [kg/kg]

   real(r8) esat
   real(r8) qsat
   real(r8) dqsdT

   real(r8) esat_b(pcols)                                 
   real(r8) qsat_b(pcols)                                 
   real(r8) dqsdT_b(pcols)                                 

   real(r8) QQmax,QQmin,QQwmin,QQimin                      ! For limiting QQ
   real(r8) cone                                           ! Number close to but smaller than 1
   real(r8) qsmall                                         ! Smallest mixing ratio considered in the macrophysics

   cone            = 0.999_r8
   qsmall          = 1.e-18_r8
   zeros(:ncol,:)  = 0._r8

   ! ------------------------------------ !
   ! Global initialization of main output !
   ! ------------------------------------ !

     s_tendout(:ncol,:)     = 0._r8
     qv_tendout(:ncol,:)    = 0._r8
     ql_tendout(:ncol,:)    = 0._r8
     qi_tendout(:ncol,:)    = 0._r8
     nl_tendout(:ncol,:)    = 0._r8
     ni_tendout(:ncol,:)    = 0._r8

     qme(:ncol,:)           = 0._r8

     cld(:ncol,:)           = 0._r8
     al_st_star(:ncol,:)    = 0._r8
     ai_st_star(:ncol,:)    = 0._r8
     ql_st_star(:ncol,:)    = 0._r8
     qi_st_star(:ncol,:)    = 0._r8

   ! --------------------------------------- !
   ! Initialization of internal 2D variables !
   ! --------------------------------------- !

     T(:ncol,:)             = 0._r8
     T1(:ncol,:)            = 0._r8
     T_0(:ncol,:)           = 0._r8
     T_05(:ncol,:)          = 0._r8
     T_prime0(:ncol,:)      = 0._r8
     T_dprime(:ncol,:)      = 0._r8
     T_star(:ncol,:)        = 0._r8

     qv(:ncol,:)            = 0._r8
     qv1(:ncol,:)           = 0._r8
     qv_0(:ncol,:)          = 0._r8
     qv_05(:ncol,:)         = 0._r8
     qv_prime0(:ncol,:)     = 0._r8
     qv_dprime(:ncol,:)     = 0._r8
     qv_star(:ncol,:)       = 0._r8

     ql(:ncol,:)            = 0._r8
     ql1(:ncol,:)           = 0._r8
     ql_0(:ncol,:)          = 0._r8
     ql_05(:ncol,:)         = 0._r8
     ql_prime0(:ncol,:)     = 0._r8
     ql_dprime(:ncol,:)     = 0._r8
     ql_star(:ncol,:)       = 0._r8

     qi(:ncol,:)            = 0._r8
     qi1(:ncol,:)           = 0._r8
     qi_0(:ncol,:)          = 0._r8
     qi_05(:ncol,:)         = 0._r8
     qi_prime0(:ncol,:)     = 0._r8
     qi_dprime(:ncol,:)     = 0._r8
     qi_star(:ncol,:)       = 0._r8

     nl(:ncol,:)            = 0._r8
     nl1(:ncol,:)           = 0._r8
     nl_0(:ncol,:)          = 0._r8
     nl_05(:ncol,:)         = 0._r8
     nl_prime0(:ncol,:)     = 0._r8
     nl_dprime(:ncol,:)     = 0._r8
     nl_star(:ncol,:)       = 0._r8

     ni(:ncol,:)            = 0._r8
     ni1(:ncol,:)           = 0._r8
     ni_0(:ncol,:)          = 0._r8
     ni_05(:ncol,:)         = 0._r8
     ni_prime0(:ncol,:)     = 0._r8
     ni_dprime(:ncol,:)     = 0._r8
     ni_star(:ncol,:)       = 0._r8

     a_st(:ncol,:)          = 0._r8
     a_st_0(:ncol,:)        = 0._r8
     a_st_star(:ncol,:)     = 0._r8

     al_st(:ncol,:)         = 0._r8
     al_st_0(:ncol,:)       = 0._r8
     al_st_nc(:ncol,:)      = 0._r8

     ai_st(:ncol,:)         = 0._r8
     ai_st_0(:ncol,:)       = 0._r8
     ai_st_nc(:ncol,:)      = 0._r8

     ql_st(:ncol,:)         = 0._r8
     ql_st_0(:ncol,:)       = 0._r8

     qi_st(:ncol,:)         = 0._r8
     qi_st_0(:ncol,:)       = 0._r8

 ! Cumulus properties 

     dacudt(:ncol,:)        = 0._r8
     a_cu(:ncol,:)          = 0._r8

 ! Adjustment tendency in association with 'positive_moisture'

     Tten_pwi1(:ncol,:)     = 0._r8
     qvten_pwi1(:ncol,:)    = 0._r8
     qlten_pwi1(:ncol,:)    = 0._r8
     qiten_pwi1(:ncol,:)    = 0._r8
     nlten_pwi1(:ncol,:)    = 0._r8
     niten_pwi1(:ncol,:)    = 0._r8

     Tten_pwi2(:ncol,:)     = 0._r8
     qvten_pwi2(:ncol,:)    = 0._r8
     qlten_pwi2(:ncol,:)    = 0._r8
     qiten_pwi2(:ncol,:)    = 0._r8
     nlten_pwi2(:ncol,:)    = 0._r8
     niten_pwi2(:ncol,:)    = 0._r8

     A_T_adj(:ncol,:)       = 0._r8
     A_qv_adj(:ncol,:)      = 0._r8
     A_ql_adj(:ncol,:)      = 0._r8
     A_qi_adj(:ncol,:)      = 0._r8
     A_nl_adj(:ncol,:)      = 0._r8
     A_ni_adj(:ncol,:)      = 0._r8

     qvadj   (:ncol,:)      = 0._r8
     qladj   (:ncol,:)      = 0._r8
     qiadj   (:ncol,:)      = 0._r8

 ! Adjustment tendency in association with 'instratus_condensate'

     QQw1(:ncol,:)          = 0._r8
     QQi1(:ncol,:)          = 0._r8
     QQw2(:ncol,:)          = 0._r8
     QQi2(:ncol,:)          = 0._r8

     QQnl1(:ncol,:)         = 0._r8
     QQni1(:ncol,:)         = 0._r8
     QQnl2(:ncol,:)         = 0._r8
     QQni2(:ncol,:)         = 0._r8

     QQnl(:ncol,:)          = 0._r8
     QQni(:ncol,:)          = 0._r8

 ! Macrophysical process tendency variables

     QQ(:ncol,:)            = 0._r8
     QQw(:ncol,:)           = 0._r8
     QQi(:ncol,:)           = 0._r8
     QQnl(:ncol,:)          = 0._r8
     QQni(:ncol,:)          = 0._r8
     ACnl(:ncol,:)          = 0._r8
     ACni(:ncol,:)          = 0._r8

     QQw_prev(:ncol,:)      = 0._r8
     QQi_prev(:ncol,:)      = 0._r8
     QQnl_prev(:ncol,:)     = 0._r8
     QQni_prev(:ncol,:)     = 0._r8

     QQw_prog(:ncol,:)      = 0._r8
     QQi_prog(:ncol,:)      = 0._r8
     QQnl_prog(:ncol,:)     = 0._r8
     QQni_prog(:ncol,:)     = 0._r8

     QQ_final(:ncol,:)      = 0._r8                        
     QQw_final(:ncol,:)     = 0._r8                  
     QQi_final(:ncol,:)     = 0._r8           
     QQn_final(:ncol,:)     = 0._r8    
     QQnl_final(:ncol,:)    = 0._r8
     QQni_final(:ncol,:)    = 0._r8

     QQ_all(:ncol,:)        = 0._r8
     QQw_all(:ncol,:)       = 0._r8
     QQi_all(:ncol,:)       = 0._r8
     QQn_all(:ncol,:)       = 0._r8
     QQnl_all(:ncol,:)      = 0._r8
     QQni_all(:ncol,:)      = 0._r8

 ! Coefficient for computing QQ and related processes

     U(:ncol,:)             = 0._r8
     U_nc(:ncol,:)          = 0._r8
     G_nc(:ncol,:)          = 0._r8
     F_nc(:ncol,:)          = 0._r8

 ! Other

     qmin1(:ncol,:)         = 0._r8
     qmin2(:ncol,:)         = 0._r8
     qmin3(:ncol,:)         = 0._r8

     esat_b(:ncol)          = 0._r8     
     qsat_b(:ncol)          = 0._r8    
     dqsdT_b(:ncol)         = 0._r8 

     esat_a(:ncol)          = 0._r8     
     qsat_a(:ncol)          = 0._r8    
     dqsdT_a(:ncol)         = 0._r8 
     Twb_aw(:ncol,:pver)    = 0._r8
     qvwb_aw(:ncol,:pver)   = 0._r8

   ! ---------------- !
   ! Main computation ! 
   ! ---------------- !

   ! ---------------------------------- !
   ! Compute cumulus-related properties ! 
   ! ---------------------------------- !

   dacudt(:ncol,:)  =  (a_cu0(:ncol,:) - a_cud(:ncol,:))/dt

   ! ---------------------------------------------------------------------- !
   ! Check if input non-cumulus pixels satisfie a non-negative constraint.  !
   ! If not, force all water vapor substances to be positive in all layers. !
   ! We should use 'old' cumulus properties for this routine.               !                
   ! ---------------------------------------------------------------------- !

   T1(:ncol,:)    =  T0(:ncol,:) 
   qv1(:ncol,:)   = qv0(:ncol,:) 
   ql1(:ncol,:)   = ql0(:ncol,:) 
   qi1(:ncol,:)   = qi0(:ncol,:) 
   nl1(:ncol,:)   = nl0(:ncol,:) 
   ni1(:ncol,:)   = ni0(:ncol,:) 

   qmin1(:ncol,:) = qmin(1)
   qmin2(:ncol,:) = qmin(2)
   qmin3(:ncol,:) = qmin(3)

   call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, & 
                           qv1, ql1, qi1, T1, qvten_pwi1, qlten_pwi1, qiten_pwi1, Tten_pwi1 )

   do k = 1, pver
   do i = 1, ncol
      if( ql1(i,k) .lt. qsmall ) then
          nlten_pwi1(i,k) = -nl1(i,k)/dt
          nl1(i,k)        = 0._r8
      endif 
      if( qi1(i,k) .lt. qsmall ) then
          niten_pwi1(i,k) = -ni1(i,k)/dt
          ni1(i,k)        = 0._r8
      endif 
   enddo
   enddo

   ! ------------------------------------------------------------- !
   ! Impose 'in-stratus condensate amount constraint'              !
   ! such that it is bounded by two limiting values.               !      
   ! This should also use 'old' cumulus properties since it is     !
   ! before applying external forcings.                            ! 
   ! Below 'QQw1,QQi1' are effective adjustive condensation        ! 
   ! Although this process also involves freezing of cloud         !
   ! liquid into ice, they can be and only can be expressed        !
   ! in terms of effective condensation.                           !
   ! ------------------------------------------------------------- !

   do k = 1, pver
      call instratus_condensate( lchnk, ncol, k,                                   &
                                 p(:,k), T1(:,k), qv1(:,k), ql1(:,k), qi1(:,k),    &
                                 a_cud(:,k), zeros(:,k), zeros(:,k),               &
                                 zeros(:,k), zeros(:,k), zeros(:,k),               &
                                 landfrac, snowh,                                  &
                                 T_0(:,k), qv_0(:,k), ql_0(:,k), qi_0(:,k),        & 
                                 al_st_0(:,k), ai_st_0(:,k), ql_st_0(:,k), qi_st_0(:,k) )
      a_st_0(:ncol,k) = max(al_st_0(:ncol,k),ai_st_0(:ncol,k))
      QQw1(:ncol,k)   = (ql_0(:ncol,k) - ql1(:ncol,k))/dt
      QQi1(:ncol,k)   = (qi_0(:ncol,k) - qi1(:ncol,k))/dt
      ! -------------------------------------------------- !
      ! Reduce droplet concentration if evaporation occurs !
      ! Set a limit such that negative state not happens.  ! 
      ! -------------------------------------------------- !
      do i = 1, ncol
         if( QQw1(i,k) .le. 0._r8 ) then
             if( ql1(i,k) .gt. qsmall ) then
                 QQnl1(i,k) = QQw1(i,k)*nl1(i,k)/ql1(i,k)
                 QQnl1(i,k) = min(0._r8,cone*max(QQnl1(i,k),-nl1(i,k)/dt))
             else
                 QQnl1(i,k) = 0._r8
             endif  
         endif 
         if( QQi1(i,k) .le. 0._r8 ) then
             if( qi1(i,k) .gt. qsmall ) then
                 QQni1(i,k) = QQi1(i,k)*ni1(i,k)/qi1(i,k)
                 QQni1(i,k) = min(0._r8,cone*max(QQni1(i,k),-ni1(i,k)/dt))
             else
                 QQni1(i,k) = 0._r8
             endif  
         endif 
      enddo
   enddo
   nl_0(:ncol,:) = max(0._r8,nl1(:ncol,:)+QQnl1(:ncol,:)*dt) 
   ni_0(:ncol,:) = max(0._r8,ni1(:ncol,:)+QQni1(:ncol,:)*dt)

   ! ----------------------------------------------------------------------------- !
   ! Check if non-cumulus pixels of '_05' state satisfies non-negative constraint. !
   ! If not, force all water substances of '_05' state to be positive by imposing  !
   ! adjustive advection. We should use 'new' cumulus properties for this routine. !                
   ! ----------------------------------------------------------------------------- !

   T_05(:ncol,:)  =  T_0(:ncol,:) + (  A_T(:ncol,:) +  C_T(:ncol,:) ) * dt
   qv_05(:ncol,:) = qv_0(:ncol,:) + ( A_qv(:ncol,:) + C_qv(:ncol,:) ) * dt
   ql_05(:ncol,:) = ql_0(:ncol,:) + ( A_ql(:ncol,:) + C_ql(:ncol,:) ) * dt
   qi_05(:ncol,:) = qi_0(:ncol,:) + ( A_qi(:ncol,:) + C_qi(:ncol,:) ) * dt 
   nl_05(:ncol,:) = max(0._r8, nl_0(:ncol,:) + ( A_nl(:ncol,:) + C_nl(:ncol,:) ) * dt )
   ni_05(:ncol,:) = max(0._r8, ni_0(:ncol,:) + ( A_ni(:ncol,:) + C_ni(:ncol,:) ) * dt )

   call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp, & 
                           qv_05, ql_05, qi_05, T_05, A_qv_adj, A_ql_adj, A_qi_adj, A_T_adj )

   ! -------------------------------------------------------------- !
   ! Define reference state at the first iteration. This will be    !
   ! continuously updated within the iteration loop below.          !
   ! While equlibrium state properties are already output from the  !
   ! 'instratus_condensate', they will be re-computed within the    !
   ! each iteration process. At the first iteration, they will      !
   ! produce exactly identical results. Note that except at the     !
   ! very first iteration iter = 1, we must use updated cumulus     !
   ! properties at all the other iteration processes. Even at the   !
   ! first iteration, we should use updated cumulus properties      !
   ! when computing limiters for (Q,P,E).                           !
   ! -------------------------------------------------------------- !

   ! -------------------------------------------------------------- !
   ! Define variables at the reference state of the first iteration !
   ! -------------------------------------------------------------- !

   T(:ncol,:)     = T_0(:ncol,:)
   qv(:ncol,:)    = qv_0(:ncol,:)
   ql(:ncol,:)    = ql_0(:ncol,:)
   qi(:ncol,:)    = qi_0(:ncol,:)
   al_st(:ncol,:) = al_st_0(:ncol,:)
   ai_st(:ncol,:) = ai_st_0(:ncol,:)
   a_st(:ncol,:)  = a_st_0(:ncol,:)
   ql_st(:ncol,:) = ql_st_0(:ncol,:)
   qi_st(:ncol,:) = qi_st_0(:ncol,:)
   nl(:ncol,:)    = nl_0(:ncol,:)
   ni(:ncol,:)    = ni_0(:ncol,:)

   ! -------------------------- !
   ! Main iterative computation !
   ! -------------------------- !

   do iter = 1, niter

      ! ------------------------------------------ !
      ! Initialize array within the iteration loop !
      ! ------------------------------------------ !

      QQ(:,:)         = 0._r8
      QQw(:,:)        = 0._r8
      QQi(:,:)        = 0._r8
      QQnl(:,:)       = 0._r8
      QQni(:,:)       = 0._r8 
      QQw2(:,:)       = 0._r8
      QQi2(:,:)       = 0._r8
      QQnl2(:,:)      = 0._r8
      QQni2(:,:)      = 0._r8
      nlten_pwi2(:,:) = 0._r8
      niten_pwi2(:,:) = 0._r8
      ACnl(:,:)       = 0._r8
      ACni(:,:)       = 0._r8 
      aa(:,:)         = 0._r8
      bb(:,:)         = 0._r8

      call findsp_water(lchnk,ncol,qv_05(:,:),T_05(:,:),p(:,:),Twb_aw(:,:),qvwb_aw(:,:))

      do k = 1, pver

      call vqsatd2_water(T_05(1:ncol,k),p(1:ncol,k),esat_a(1:ncol),qsat_a(1:ncol),dqsdT_a(1:ncol),ncol)
      call vqsatd2_water(T(1:ncol,k),p(1:ncol,k),esat_b(1:ncol),qsat_b(1:ncol),dqsdT_b(1:ncol),ncol)

      if( iter .eq. 1 ) then
          a_cu(:ncol,k) = a_cud(:ncol,k)
      else
          a_cu(:ncol,k) = a_cu0(:ncol,k)
      endif
      do i = 1, ncol
         U(i,k)    =  qv(i,k)/qsat_b(i)
         U_nc(i,k) =  U(i,k)
      enddo
      if( CAMstfrac ) then
          call astG_RHU(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol)
      else
          call astG_PDF(U_nc(:,k),p(:,k),qv(:,k),landfrac(:),snowh(:),al_st_nc(:,k),G_nc(:,k),ncol)
      endif
      call aist_vector(qv(:,k),T(:,k),p(:,k),qi(:,k),landfrac(:),snowh(:),ai_st_nc(:,k),ncol)
      ai_st(:ncol,k)  =  (1._r8-a_cu(:ncol,k))*ai_st_nc(:ncol,k)
      al_st(:ncol,k)  =  (1._r8-a_cu(:ncol,k))*al_st_nc(:ncol,k)
      a_st(:ncol,k)   =  max(al_st(:ncol,k),ai_st(:ncol,k))  

      do i = 1, ncol

         ! -------------------------------------------------------- !
         ! Compute basic thermodynamic coefficients for computing Q !
         ! -------------------------------------------------------- !

         alpha  =  1._r8/qsat_b(i)
         beta   =  dqsdt_b(i)*(qv(i,k)/qsat_b(i)**2)
         betast =  alpha*dqsdT_b(i) 
         gammal =  alpha + (hlatv/cp)*beta
         gammai =  alpha + ((hlatv+hlatf)/cp)*beta
         gammaQ =  alpha + (hlatv/cp)*beta
         deltal =  1._r8 + a_st(i,k)*(hlatv/cp)*(betast/alpha)
         deltai =  1._r8 + a_st(i,k)*((hlatv+hlatf)/cp)*(betast/alpha)
         A_Tc   =  A_T(i,k)+A_T_adj(i,k)-(hlatv/cp)*(A_ql(i,k)+A_ql_adj(i,k))-((hlatv+hlatf)/cp)*(A_qi(i,k)+A_qi_adj(i,k))
         A_qt   =  A_qv(i,k) + A_qv_adj(i,k) + A_ql(i,k) + A_ql_adj(i,k) + A_qi(i,k) + A_qi_adj(i,k)
         C_Tc   =  C_T(i,k) - (hlatv/cp)*C_ql(i,k) - ((hlatv+hlatf)/cp)*C_qi(i,k)
         C_qt   =  C_qv(i,k) + C_ql(i,k) + C_qi(i,k)
         dTcdt  =  A_Tc + C_Tc
         dqtdt  =  A_qt + C_qt
       ! dqtstldt = A_qt + C_ql(i,k)/max(1.e-2_r8,al_st(i,k))                             ! Original  
       ! dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_ql(i,k)/max(1.e-2_r8,al_st(i,k)) ! New 1 on Dec.30.2009.
         dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_qlst(i,k)                        ! New 2 on Dec.30.2009.
       ! dqtstldt = A_qt + C_qt                                                           ! Original Conservative treatment
       ! dqtstldt = A_qt - A_qi(i,k) - A_qi_adj(i,k) + C_qt - C_qi(i,k)                   ! New Conservative treatment on Dec.30.2009
         dqidt = A_qi(i,k) + A_qi_adj(i,k) + C_qi(i,k) 

         anic    = max(1.e-8_r8,(1._r8-a_cu(i,k)))
         GG      = G_nc(i,k)/anic
         aa(1,1) = gammal*al_st(i,k)
         aa(1,2) = GG + gammal*cc*ql_st(i,k)          
         aa(2,1) = alpha + (hlatv/cp)*betast*al_st(i,k)
         aa(2,2) = (hlatv/cp)*betast*cc*ql_st(i,k) 
         bb(1,1) = alpha*dqtdt - beta*dTcdt - gammai*dqidt - GG*al_st_nc(i,k)*dacudt(i,k) + F_nc(i,k) 
         bb(2,1) = alpha*dqtstldt - betast*(dTcdt + ((hlatv+hlatf)/cp)*dqidt) 
         call gaussj(aa(1:2,1:2),2,2,bb(1:2,1),1,1)
         dqlstdt = bb(1,1)
         dalstdt = bb(2,1)
         QQ(i,k) = al_st(i,k)*dqlstdt + cc*ql_st(i,k)*dalstdt - ( A_ql(i,k) + A_ql_adj(i,k) + C_ql(i,k) )

       ! ------------------------------------------------------------ !
       ! Limiter for QQ                                               !
       ! Here, 'fice' should be from the reference equilibrium state  !
       ! since QQ itself is computed from the reference state.        !
       ! From the assumption used for derivation of QQ(i), it must be !
       ! that QQw(i) = QQ(i)*(1._r8-fice(i)), QQi(i) = QQ(i)*fice(i)  !  
       ! ------------------------------------------------------------ !

         if( QQ(i,k) .ge. 0._r8 ) then
             QQmax    = (qv_05(i,k) - qmin(1))/dt ! For ghost cumulus & semi-ghost ice stratus
             QQmax    = max(0._r8,QQmax) 
             QQ(i,k)  = min(QQ(i,k),QQmax)
             QQw(i,k) = QQ(i,k)
             QQi(i,k) = 0._r8 
         else
             QQmin  = 0._r8
             if( qv_05(i,k) .lt. qsat_a(i) ) QQmin = min(0._r8,cone*(qv_05(i,k)-qvwb_aw(i,k))/dt)
             QQ(i,k)  = max(QQ(i,k),QQmin)
             QQw(i,k) = QQ(i,k)
             QQi(i,k) = 0._r8
             QQwmin   = min(0._r8,-cone*ql_05(i,k)/dt)
             QQimin   = min(0._r8,-cone*qi_05(i,k)/dt)
             QQw(i,k) = min(0._r8,max(QQw(i,k),QQwmin))
             QQi(i,k) = min(0._r8,max(QQi(i,k),QQimin))
         endif

       ! -------------------------------------------------- !
       ! Reduce droplet concentration if evaporation occurs !
       ! Note 'QQnl1,QQni1' are computed from the reference !
       ! equilibrium state but limiter is from 'nl_05'.     !
       ! -------------------------------------------------- !

         if( QQw(i,k) .lt. 0._r8 ) then
             if( ql_05(i,k) .gt. qsmall ) then
                 QQnl(i,k) = QQw(i,k)*nl_05(i,k)/ql_05(i,k)
                 QQnl(i,k) = min(0._r8,cone*max(QQnl(i,k),-nl_05(i,k)/dt))
             else
                 QQnl(i,k) = 0._r8
             endif  
         endif 

         if( QQi(i,k) .lt. 0._r8 ) then
             if( qi_05(i,k) .gt. qsmall ) then
                 QQni(i,k) = QQi(i,k)*ni_05(i,k)/qi_05(i,k)
                 QQni(i,k) = min(0._r8,cone*max(QQni(i,k),-ni_05(i,k)/dt))
             else
                 QQni(i,k) = 0._r8
             endif  
         endif 

      enddo
      enddo

    ! -------------------------------------------------------------------- !
    ! Until now, we have finished computing all necessary tendencies       ! 
    ! from the equilibrium input state (T_0).                              !
    ! If ramda = 0 : fully explicit scheme                                 !
    !    ramda = 1 : fully implicit scheme                                 !
    ! Note that 'ramda = 0.5 with niter = 2' can mimic                     !
    ! -------------------------------------------------------------------- !

      if( iter .eq. 1 ) then
          QQw_prev(:ncol,:)  = QQw(:ncol,:)       
          QQi_prev(:ncol,:)  = QQi(:ncol,:)   
          QQnl_prev(:ncol,:) = QQnl(:ncol,:)       
          QQni_prev(:ncol,:) = QQni(:ncol,:)   
      endif

      QQw_prog(:ncol,:)   = ramda*QQw(:ncol,:)   + (1._r8-ramda)*QQw_prev(:ncol,:)
      QQi_prog(:ncol,:)   = ramda*QQi(:ncol,:)   + (1._r8-ramda)*QQi_prev(:ncol,:)
      QQnl_prog(:ncol,:)  = ramda*QQnl(:ncol,:)  + (1._r8-ramda)*QQnl_prev(:ncol,:)
      QQni_prog(:ncol,:)  = ramda*QQni(:ncol,:)  + (1._r8-ramda)*QQni_prev(:ncol,:)

      QQw_prev(:ncol,:)   = QQw_prog(:ncol,:)
      QQi_prev(:ncol,:)   = QQi_prog(:ncol,:)
      QQnl_prev(:ncol,:)  = QQnl_prog(:ncol,:)
      QQni_prev(:ncol,:)  = QQni_prog(:ncol,:)

    ! -------------------------------------------------------- !
    ! Compute final prognostic state on which final diagnostic !
    ! in-stratus condensate adjustment is applied in the below.!
    ! Important : I must check whether there are any external  !  
    !             advective forcings of 'A_nl(i,k),A_ni(i,k)'. !
    !             Even they are (i.e., advection of aerosol),  !
    !             actual droplet activation will be performd   !
    !             in microphysics, so it will be completely    !
    !             reasonable to 'A_nl(i,k)=A_ni(i,k)=0'.       !
    ! -------------------------------------------------------- !

    do k = 1, pver
    do i = 1, ncol
       T_prime0(i,k)  = T_0(i,k)  + dt*( A_T(i,k)  +  A_T_adj(i,k) +  C_T(i,k) + (hlatv*QQw_prog(i,k)+(hlatv+hlatf)*QQi_prog(i,k))/cp )
       qv_prime0(i,k) = qv_0(i,k) + dt*( A_qv(i,k) + A_qv_adj(i,k) + C_qv(i,k) - QQw_prog(i,k) - QQi_prog(i,k) )
       ql_prime0(i,k) = ql_0(i,k) + dt*( A_ql(i,k) + A_ql_adj(i,k) + C_ql(i,k) + QQw_prog(i,k) )
       qi_prime0(i,k) = qi_0(i,k) + dt*( A_qi(i,k) + A_qi_adj(i,k) + C_qi(i,k) + QQi_prog(i,k) )
       nl_prime0(i,k) = max(0._r8,nl_0(i,k) + dt*( A_nl(i,k) + C_nl(i,k) + QQnl_prog(i,k) ))
       ni_prime0(i,k) = max(0._r8,ni_0(i,k) + dt*( A_ni(i,k) + C_ni(i,k) + QQni_prog(i,k) ))
       if( ql_prime0(i,k) .lt. qsmall ) nl_prime0(i,k) = 0._r8
       if( qi_prime0(i,k) .lt. qsmall ) ni_prime0(i,k) = 0._r8
    enddo
    enddo

   ! -------------------------------------------------- !
   ! Perform diagnostic 'positive_moisture' constraint. !
   ! -------------------------------------------------- !

   T_dprime(:ncol,:)  =  T_prime0(:ncol,:) 
   qv_dprime(:ncol,:) = qv_prime0(:ncol,:) 
   ql_dprime(:ncol,:) = ql_prime0(:ncol,:) 
   qi_dprime(:ncol,:) = qi_prime0(:ncol,:) 
   nl_dprime(:ncol,:) = nl_prime0(:ncol,:) 
   ni_dprime(:ncol,:) = ni_prime0(:ncol,:) 

   call positive_moisture( ncol, dt, qmin1, qmin2, qmin3, dp,          & 
                           qv_dprime, ql_dprime, qi_dprime, T_dprime,  &
                           qvten_pwi2, qlten_pwi2, qiten_pwi2, Tten_pwi2 )

   do k = 1, pver
   do i = 1, ncol
      if( ql_dprime(i,k) .lt. qsmall ) then
          nlten_pwi2(i,k) = -nl_dprime(i,k)/dt
          nl_dprime(i,k)   = 0._r8
      endif 
      if( qi_dprime(i,k) .lt. qsmall ) then
          niten_pwi2(i,k) = -ni_dprime(i,k)/dt
          ni_dprime(i,k)   = 0._r8
      endif 
   enddo
   enddo

   ! -------------------------------------------------------------- !
   ! Add tendency associated with detrainment of cumulus condensate !
   ! This tendency is not used in computing Q                       !
   ! Since D_ql,D_qi,D_nl,D_ni > 0, don't need to worry about       !
   ! negative scalar.                                               !
   ! This tendency is not reflected into Fzs2, which is OK.         !
   ! -------------------------------------------------------------- !

   T_dprime(:ncol,:)   =  T_dprime(:ncol,:)  + D_T(:ncol,:) * dt 
   qv_dprime(:ncol,:)  = qv_dprime(:ncol,:) + D_qv(:ncol,:) * dt 
   ql_dprime(:ncol,:)  = ql_dprime(:ncol,:) + D_ql(:ncol,:) * dt
   qi_dprime(:ncol,:)  = qi_dprime(:ncol,:) + D_qi(:ncol,:) * dt
   nl_dprime(:ncol,:)  = nl_dprime(:ncol,:) + D_nl(:ncol,:) * dt 
   ni_dprime(:ncol,:)  = ni_dprime(:ncol,:) + D_ni(:ncol,:) * dt

   ! ---------------------------------------------------------- !
   ! Impose diagnostic upper and lower limits on the in-stratus !
   ! condensate amount. This produces a final equilibrium state !
   ! at the end of each iterative process.                      !
   ! ---------------------------------------------------------- !

   do k = 1, pver
      call instratus_condensate( lchnk          , ncol           , k              , p(:,k)        , &
                                 T_dprime(:,k)  , qv_dprime(:,k) , ql_dprime(:,k) , qi_dprime(:,k), &
                                 a_cu0(:,k)     , zeros(:,k)     , zeros(:,k)     ,                 & 
                                 zeros(:,k)     , zeros(:,k)     , zeros(:,k)     ,                 &
                                 landfrac       , snowh          ,                                  &
                                 T_star(:,k)    , qv_star(:,k)   , ql_star(:,k)   , qi_star(:,k)  , & 
                                 al_st_star(:,k), ai_st_star(:,k), ql_st_star(:,k), qi_st_star(:,k) )
      a_st_star(:ncol,k)  = max(al_st_star(:ncol,k),ai_st_star(:ncol,k))
      QQw2(:ncol,k) = (ql_star(:ncol,k) - ql_dprime(:ncol,k))/dt
      QQi2(:ncol,k) = (qi_star(:ncol,k) - qi_dprime(:ncol,k))/dt
      ! -------------------------------------------------- !
      ! Reduce droplet concentration if evaporation occurs !
      ! -------------------------------------------------- !
      do i = 1, ncol
         if( QQw2(i,k) .le. 0._r8 ) then
             if( ql_dprime(i,k) .ge. qsmall ) then
                 QQnl2(i,k) = QQw2(i,k)*nl_dprime(i,k)/ql_dprime(i,k)
                 QQnl2(i,k) = min(0._r8,cone*max(QQnl2(i,k),-nl_dprime(i,k)/dt))
             else
                 QQnl2(i,k) = 0._r8
             endif  
         endif 
         if( QQi2(i,k) .le. 0._r8 ) then
             if( qi_dprime(i,k) .gt. qsmall ) then
                 QQni2(i,k) = QQi2(i,k)*ni_dprime(i,k)/qi_dprime(i,k)
                 QQni2(i,k) = min(0._r8,cone*max(QQni2(i,k),-ni_dprime(i,k)/dt))
             else
                 QQni2(i,k) = 0._r8
             endif  
         endif 
      enddo
   enddo
   nl_star(:ncol,:) = max(0._r8,nl_dprime(:ncol,:)+QQnl2(:ncol,:)*dt) 
   ni_star(:ncol,:) = max(0._r8,ni_dprime(:ncol,:)+QQni2(:ncol,:)*dt)

   ! ------------------------------------------ !
   ! Final adjustment of droplet concentration. !
   ! Set # to zero if there is no cloud.        !
   ! ------------------------------------------ !

   do k = 1, pver
   do i = 1, ncol 
      if( ql_star(i,k) .lt. qsmall ) then
          ACnl(i,k) = - nl_star(i,k)/dt
          nl_star(i,k) = 0._r8
      endif
      if( qi_star(i,k) .lt. qsmall ) then
          ACni(i,k) = - ni_star(i,k)/dt
          ni_star(i,k) = 0._r8
      endif
   enddo
   enddo

   ! ----------------------------------------------------- !
   ! Define equilibrium reference state for next iteration !
   ! ----------------------------------------------------- !

   T(:ncol,:)     = T_star(:ncol,:)
   qv(:ncol,:)    = qv_star(:ncol,:)
   ql(:ncol,:)    = ql_star(:ncol,:)
   qi(:ncol,:)    = qi_star(:ncol,:)
   al_st(:ncol,:) = al_st_star(:ncol,:)
   ai_st(:ncol,:) = ai_st_star(:ncol,:)
   a_st(:ncol,:)  = a_st_star(:ncol,:)
   ql_st(:ncol,:) = ql_st_star(:ncol,:)
   qi_st(:ncol,:) = qi_st_star(:ncol,:)
   nl(:ncol,:)    = nl_star(:ncol,:)
   ni(:ncol,:)    = ni_star(:ncol,:)

   enddo ! End of 'iter' prognostic iterative computation

   ! ------------------------------------------------------------------------ !
   ! Compute final tendencies of main output variables and diagnostic outputs !
   ! Note that the very input state [T0,qv0,ql0,qi0] are                      !
   ! marched to [T_star,qv_star,ql_star,qi_star] with equilibrium             !
   ! stratus informations of [a_st_star,ql_st_star,qi_st_star] by             !
   ! below final tendencies and [A_T,A_qv,A_ql,A_qi]                          !
   ! ------------------------------------------------------------------------ !

   ! ------------------ !
   ! Process tendencies !
   ! ------------------ !

   QQw_final(:ncol,:)  = QQw_prog(:ncol,:)
   QQi_final(:ncol,:)  = QQi_prog(:ncol,:)
   QQ_final(:ncol,:)   = QQw_final(:ncol,:) + QQi_final(:ncol,:)
   QQw_all(:ncol,:)    = QQw_prog(:ncol,:)  + QQw1(:ncol,:) + QQw2(:ncol,:) + qlten_pwi1(:ncol,:) + qlten_pwi2(:ncol,:) + A_ql_adj(:ncol,:)
   QQi_all(:ncol,:)    = QQi_prog(:ncol,:)  + QQi1(:ncol,:) + QQi2(:ncol,:) + qiten_pwi1(:ncol,:) + qiten_pwi2(:ncol,:) + A_qi_adj(:ncol,:)
   QQ_all(:ncol,:)     = QQw_all(:ncol,:)   + QQi_all(:ncol,:)
   QQnl_final(:ncol,:) = QQnl_prog(:ncol,:)
   QQni_final(:ncol,:) = QQni_prog(:ncol,:)
   QQn_final(:ncol,:)  = QQnl_final(:ncol,:) + QQni_final(:ncol,:)
   QQnl_all(:ncol,:)   = QQnl_prog(:ncol,:)  + QQnl1(:ncol,:) + QQnl2(:ncol,:) + nlten_pwi1(:ncol,:) + nlten_pwi2(:ncol,:) + ACnl(:ncol,:) + A_nl_adj(:ncol,:)
   QQni_all(:ncol,:)   = QQni_prog(:ncol,:)  + QQni1(:ncol,:) + QQni2(:ncol,:) + niten_pwi1(:ncol,:) + niten_pwi2(:ncol,:) + ACni(:ncol,:) + A_ni_adj(:ncol,:)
   QQn_all(:ncol,:)    = QQnl_all(:ncol,:)   + QQni_all(:ncol,:)
   qme(:ncol,:)        = QQ_final(:ncol,:)   
   qvadj(:ncol,:)      = qvten_pwi1(:ncol,:) + qvten_pwi2(:ncol,:) + A_qv_adj(:ncol,:)
   qladj(:ncol,:)      = qlten_pwi1(:ncol,:) + qlten_pwi2(:ncol,:) + A_ql_adj(:ncol,:)
   qiadj(:ncol,:)      = qiten_pwi1(:ncol,:) + qiten_pwi2(:ncol,:) + A_qi_adj(:ncol,:)
   qllim(:ncol,:)      = QQw1      (:ncol,:) + QQw2      (:ncol,:)
   qilim(:ncol,:)      = QQi1      (:ncol,:) + QQi2      (:ncol,:)

   ! ----------------- !
   ! Output tendencies !
   ! ----------------- !

   s_tendout(:ncol,:)  = cp*( T_star(:ncol,:)  -  T0(:ncol,:) )/dt - cp*(A_T(:ncol,:)+C_T(:ncol,:))
   qv_tendout(:ncol,:) =    ( qv_star(:ncol,:) - qv0(:ncol,:) )/dt - (A_qv(:ncol,:)+C_qv(:ncol,:))
   ql_tendout(:ncol,:) =    ( ql_star(:ncol,:) - ql0(:ncol,:) )/dt - (A_ql(:ncol,:)+C_ql(:ncol,:))
   qi_tendout(:ncol,:) =    ( qi_star(:ncol,:) - qi0(:ncol,:) )/dt - (A_qi(:ncol,:)+C_qi(:ncol,:))
   nl_tendout(:ncol,:) =    ( nl_star(:ncol,:) - nl0(:ncol,:) )/dt - (A_nl(:ncol,:)+C_nl(:ncol,:))
   ni_tendout(:ncol,:) =    ( ni_star(:ncol,:) - ni0(:ncol,:) )/dt - (A_ni(:ncol,:)+C_ni(:ncol,:))

   ! ------------------ !
   ! Net cloud fraction !
   ! ------------------ !

   cld(:ncol,:) = a_st_star(:ncol,:) + a_cu0(:ncol,:)

   ! --------------------------------- !
   ! Updated grid-mean state variables !
   ! --------------------------------- !

   T0(:ncol,:)  = T_star(:ncol,:)
   qv0(:ncol,:) = qv_star(:ncol,:)
   ql0(:ncol,:) = ql_star(:ncol,:)
   qi0(:ncol,:) = qi_star(:ncol,:)
   nl0(:ncol,:) = nl_star(:ncol,:)
   ni0(:ncol,:) = ni_star(:ncol,:)

   return
   end subroutine mmacro_pcond

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine instratus_condensate( lchnk, ncol, k,                      &  
                                    p_in, T0_in, qv0_in, ql0_in, qi0_in, & 
                                    a_dc_in, ql_dc_in, qi_dc_in,         &
                                    a_sc_in, ql_sc_in, qi_sc_in,         & 
                                    landfrac, snowh,                     &
                                    T_out, qv_out, ql_out, qi_out,       &
                                    al_st_out, ai_st_out, ql_st_out, qi_st_out )

   ! ------------------------------------------------------- !
   ! Diagnostically force in-stratus condensate to be        ! 
   ! in the range of 'qlst_min < qc_st < qlst_max'           !
   ! whenever stratus exists in the equilibrium state        !
   ! ------------------------------------------------------- !

   use time_manager,  only: is_first_step, get_nstep

   implicit none

   integer,  intent(in)  :: lchnk                ! Chunk identifier
   integer,  intent(in)  :: ncol                 ! Number of atmospheric columns
   integer,  intent(in)  :: k                    ! Layer index

   real(r8), intent(in)  :: p_in(pcols)          ! Pressure [Pa]
   real(r8), intent(in)  :: T0_in(pcols)         ! Temperature [K]
   real(r8), intent(in)  :: qv0_in(pcols)        ! Grid-mean water vapor [kg/kg]
   real(r8), intent(in)  :: ql0_in(pcols)        ! Grid-mean LWC [kg/kg]
   real(r8), intent(in)  :: qi0_in(pcols)        ! Grid-mean IWC [kg/kg]

   real(r8), intent(in)  :: a_dc_in(pcols)       ! Deep cumulus cloud fraction
   real(r8), intent(in)  :: ql_dc_in(pcols)      ! In-deep cumulus LWC [kg/kg]
   real(r8), intent(in)  :: qi_dc_in(pcols)      ! In-deep cumulus IWC [kg/kg]
   real(r8), intent(in)  :: a_sc_in(pcols)       ! Shallow cumulus cloud fraction
   real(r8), intent(in)  :: ql_sc_in(pcols)      ! In-shallow cumulus LWC [kg/kg]
   real(r8), intent(in)  :: qi_sc_in(pcols)      ! In-shallow cumulus IWC [kg/kg]

   real(r8), intent(in)  :: landfrac(pcols)      ! Land fraction
   real(r8), intent(in)  :: snowh(pcols)         ! Snow depth (liquid water equivalent)

   real(r8), intent(out) :: T_out(pcols)         ! Temperature [K]
   real(r8), intent(out) :: qv_out(pcols)        ! Grid-mean water vapor [kg/kg]
   real(r8), intent(out) :: ql_out(pcols)        ! Grid-mean LWC [kg/kg]
   real(r8), intent(out) :: qi_out(pcols)        ! Grid-mean IWC [kg/kg]

   real(r8), intent(out) :: al_st_out(pcols)     ! Liquid stratus fraction
   real(r8), intent(out) :: ai_st_out(pcols)     ! Ice stratus fraction
   real(r8), intent(out) :: ql_st_out(pcols)     ! In-stratus LWC [kg/kg]
   real(r8), intent(out) :: qi_st_out(pcols)     ! In-stratus IWC [kg/kg]

   ! Local variables

   integer i                                     ! Column    index

   real(r8) p    
   real(r8) T0   
   real(r8) qv0    
   real(r8) ql0    
   real(r8) qi0    
   real(r8) a_dc   
   real(r8) ql_dc  
   real(r8) qi_dc  
   real(r8) a_sc   
   real(r8) ql_sc  
   real(r8) qi_sc  
   real(r8) esat0  
   real(r8) qsat0  
   real(r8) dqsdT0 
   real(r8) U0     
   real(r8) U0_nc  
   real(r8) G0_nc
   real(r8) al0_st_nc            
   real(r8) al0_st
   real(r8) ai0_st_nc            
   real(r8) ai0_st               
   real(r8) a0_st               
   real(r8) ql0_nc
   real(r8) qi0_nc
   real(r8) qc0_nc
   real(r8) ql0_st
   real(r8) qi0_st
   real(r8) qc0_st
   real(r8) T   
   real(r8) qv    
   real(r8) ql    
   real(r8) qi
   real(r8) ql_st
   real(r8) qi_st
   real(r8) esat  
   real(r8) qsat  
   real(r8) dqsdT 
   real(r8) esat_in(pcols)  
   real(r8) qsat_in(pcols)  
   real(r8) dqsdT_in(pcols) 
   real(r8) U0_in(pcols)  
   real(r8) al0_st_nc_in(pcols)
   real(r8) ai0_st_nc_in(pcols)
   real(r8) G0_nc_in(pcols)
   integer  idxmod 
   real(r8) U
   real(r8) U_nc
   real(r8) al_st_nc
   real(r8) ai_st_nc
   real(r8) G_nc
   real(r8) a_st
   real(r8) al_st
   real(r8) ai_st
   real(r8) Tmin0
   real(r8) Tmax0
   real(r8) Tmin
   real(r8) Tmax
   integer caseid

   ! ---------------- !
   ! Main Computation ! 
   ! ---------------- !

   call vqsatd2_water(T0_in(1:ncol),p_in(1:ncol),esat_in(1:ncol),qsat_in(1:ncol),dqsdT_in(1:ncol),ncol)
   U0_in(:ncol) = qv0_in(:ncol)/qsat_in(:ncol)
   if( CAMstfrac ) then
       call astG_RHU(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol)
   else
       call astG_PDF(U0_in(:),p_in(:),qv0_in(:),landfrac(:),snowh(:),al0_st_nc_in(:),G0_nc_in(:),ncol)
   endif
   call aist_vector(qv0_in(:),T0_in(:),p_in(:),qi0_in(:),landfrac(:),snowh(:),ai0_st_nc_in(:),ncol)

   do i = 1, ncol

      ! ---------------------- !
      ! Define local variables !
      ! ---------------------- !

      p   = p_in(i)

      T0  = T0_in(i)
      qv0 = qv0_in(i)
      ql0 = ql0_in(i)
      qi0 = qi0_in(i)

      a_dc  = a_dc_in(i)
      ql_dc = ql_dc_in(i)
      qi_dc = qi_dc_in(i)

      a_sc  = a_sc_in(i)
      ql_sc = ql_sc_in(i)
      qi_sc = qi_sc_in(i)

      ql_dc = 0._r8
      qi_dc = 0._r8
      ql_sc = 0._r8
      qi_sc = 0._r8

      esat  = esat_in(i) 
      qsat  = qsat_in(i) 
      dqsdT = dqsdT_in(i) 

      idxmod = 0
      caseid = -1

      ! ------------------------------------------------------------ !
      ! Force the grid-mean RH to be smaller than 1 if oversaturated !
      ! In order to be compatible with reduced 3x3 QQ, condensation  !
      ! should occur only into the liquid in gridmean_RH.            !
      ! ------------------------------------------------------------ !

      if( qv0 .gt. qsat ) then
          call gridmean_RH( lchnk, i, k, p, T0, qv0, ql0, qi0,      &
                            a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, &
                            landfrac(i), snowh(i) )
          call vqsatd2_water_single(T0,p,esat0,qsat0,dqsdT0)
          U0      = (qv0/qsat0)
          U0_nc   =  U0 
          if( CAMstfrac ) then
              call astG_RHU_single(U0_nc,p,qv0,landfrac(i),snowh(i),al0_st_nc,G0_nc)
          else
              call astG_PDF_single(U0_nc,p,qv0,landfrac(i),snowh(i),al0_st_nc,G0_nc)
          endif
          call aist_single(qv0,T0,p,qi0,landfrac(i),snowh(i),ai0_st_nc)
          ai0_st  = (1._r8-a_dc-a_sc)*ai0_st_nc
          al0_st  = (1._r8-a_dc-a_sc)*al0_st_nc
          a0_st   = max(ai0_st,al0_st)         
          idxmod  = 1 
      else
          ai0_st  = (1._r8-a_dc-a_sc)*ai0_st_nc_in(i)
          al0_st  = (1._r8-a_dc-a_sc)*al0_st_nc_in(i)
      endif    
      a0_st   = max(ai0_st,al0_st)         

      ! ----------------------- ! 
      ! Handling of input state !
      ! ----------------------- !

      ql0_nc  = max(0._r8,ql0-a_dc*ql_dc-a_sc*ql_sc)
      qi0_nc  = max(0._r8,qi0-a_dc*qi_dc-a_sc*qi_sc)
      qc0_nc  = ql0_nc + qi0_nc 

      Tmin0 = T0 - (hlatv/cp)*ql0
      Tmax0 = T0 + ((hlatv+hlatf)/cp)*qv0

      ! ------------------------------------------------------------- !
      ! Do nothing and just exit if generalized in-stratus condensate !
      ! condition is satisfied. This includes the case I.             !
      ! For 4x4 liquid stratus, a0_st --> al0_st.                     ! 
      ! ------------------------------------------------------------- !
      if( ( ql0_nc .ge. qlst_min*al0_st ) .and. ( ql0_nc .le. qlst_max*al0_st ) ) then

          ! ------------------ !
          ! This is the case I !
          ! ------------------ ! 
             T = T0
             qv = qv0
             ql = ql0
             qi = qi0
             caseid = 0
             goto 10
      else
         ! ----------------------------- !
         ! This is case II : Dense Cloud !
         ! ----------------------------- !   
         if( al0_st .eq. 0._r8 .and. ql0_nc .gt. 0._r8 ) then
             ! ------------------------------------- !
             ! Compute hypothetical full evaporation !
             ! ------------------------------------- !
             T  = Tmin0
             qv = qv0 + ql0 
             call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
             U  = qv/qsat
             U_nc = U  
             if( CAMstfrac ) then
                 call astG_RHU_single(U_nc,p,qv,landfrac(i),snowh(i),al_st_nc,G_nc)
             else
                 call astG_PDF_single(U_nc,p,qv,landfrac(i),snowh(i),al_st_nc,G_nc)
             endif
             al_st = (1._r8-a_dc-a_sc)*al_st_nc  
             caseid = 0

             if( al_st .eq. 0._r8 ) then
                 ql = 0._r8
                 qi = qi0
                 idxmod = 1
                 caseid = 1
                 goto 10
             else
                 ! ------------------------------------------- !
                 ! Evaporate until qc_st decreases to qlst_max !
                 ! ------------------------------------------- !
                 Tmin = Tmin0
                 Tmax = T0 
                 call instratus_core( lchnk, i, k, p,                              &
                                      T0, qv0, ql0, 0._r8,                         &
                                      a_dc, ql_dc, qi_dc,                          &
                                      a_sc, ql_sc, qi_sc, ai0_st,                  &
                                      qlst_max, Tmin, Tmax, landfrac(i), snowh(i), &
                                      T, qv, ql, qi )   
                 idxmod = 1
                 caseid = 2
                 goto 10
             endif
         ! ------------------------------ !
         ! This is case III : Empty Cloud !
         ! ------------------------------ !  
         elseif( al0_st .gt. 0._r8 .and. ql0_nc .eq. 0._r8 ) then
              ! ------------------------------------------ ! 
              ! Condense until qc_st increases to qlst_min !
              ! ------------------------------------------ !
              Tmin = Tmin0
              Tmax = Tmax0  
              call instratus_core( lchnk, i, k, p,                              &
                                   T0, qv0, ql0, 0._r8,                         &
                                   a_dc, ql_dc, qi_dc,                          &
                                   a_sc, ql_sc, qi_sc, ai0_st,                  &
                                   qlst_min, Tmin, Tmax, landfrac(i), snowh(i), &
                                   T, qv, ql, qi )   
              idxmod = 1 
              caseid = 3
              goto 10
         ! --------------- !
         ! This is case IV !
         ! --------------- !   
         elseif( al0_st .gt. 0._r8 .and. ql0_nc .gt. 0._r8 ) then

             if( ql0_nc .gt. qlst_max*al0_st ) then
                 ! --------------------------------------- !
                 ! Evaporate until qc_st drops to qlst_max !
                 ! --------------------------------------- !
                 Tmin = Tmin0
                 Tmax = Tmax0
                 call instratus_core( lchnk, i, k, p,                              &
                                      T0, qv0, ql0, 0._r8,                         &
                                      a_dc, ql_dc, qi_dc,                          &
                                      a_sc, ql_sc, qi_sc, ai0_st,                  &
                                      qlst_max, Tmin, Tmax, landfrac(i), snowh(i), &
                                      T, qv, ql, qi )   
                 idxmod = 1
                 caseid = 4
                 goto 10
             elseif( ql0_nc .lt. qlst_min*al0_st ) then
                 ! -------------------------------------------- !
                 ! Condensate until qc_st increases to qlst_min !
                 ! -------------------------------------------- !
                 Tmin = Tmin0
                 Tmax = Tmax0 
                 call instratus_core( lchnk, i, k, p,                              &
                                      T0, qv0, ql0, 0._r8,                         &
                                      a_dc, ql_dc, qi_dc,                          &
                                      a_sc, ql_sc, qi_sc, ai0_st,                  &
                                      qlst_min, Tmin, Tmax, landfrac(i), snowh(i), & 
                                      T, qv, ql, qi )   
                 idxmod = 1
                 caseid = 5
                 goto 10
             else
                 ! ------------------------------------------------ !
                 ! This case should not happen. Issue error message !
                 ! ------------------------------------------------ !
                 write(iulog,*) 'Impossible case1 in instratus_condensate' 
                 call endrun
             endif
         ! ------------------------------------------------ !                   
         ! This case should not happen. Issue error message !
         ! ------------------------------------------------ !    
         else
             write(iulog,*) 'Impossible case2 in instratus_condensate' 
             write(iulog,*)  al0_st, a_sc, a_dc
             write(iulog,*)  1000*ql0_nc, 1000*(ql0+qi0)
             call endrun
         endif
      endif

10 continue   

   ! -------------------------------------------------- !
   ! Force final energy-moisture conserving consistency !
   ! -------------------------------------------------- !

     qi = qi0

     if( idxmod .eq. 1 ) then
         call aist_single(qv,T,p,qi,landfrac(i),snowh(i),ai_st_nc)
         ai_st = (1._r8-a_dc-a_sc)*ai_st_nc
         call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
         U     = (qv/qsat)
         U_nc  =  U
         if( CAMstfrac ) then
             call astG_RHU_single(U_nc,p,qv,landfrac(i),snowh(i),al_st_nc,G_nc)
         else
             call astG_PDF_single(U_nc,p,qv,landfrac(i),snowh(i),al_st_nc,G_nc)
         endif
         al_st = (1._r8-a_dc-a_sc)*al_st_nc
     else
         ai_st  = (1._r8-a_dc-a_sc)*ai0_st_nc_in(i)
         al_st  = (1._r8-a_dc-a_sc)*al0_st_nc_in(i)
     endif

     a_st  = max(ai_st,al_st)

     if( al_st .eq. 0._r8 ) then
         ql_st = 0._r8
     else
         ql_st = ql/al_st
         ql_st = min(qlst_max,max(qlst_min,ql_st)) ! PJR
     endif
     if( ai_st .eq. 0._r8 ) then
         qi_st = 0._r8
     else
         qi_st = qi/ai_st
     endif

     qi    = ai_st*qi_st
     ql    = al_st*ql_st

     T     = T0 - (hlatv/cp)*(ql0-ql) - ((hlatv+hlatf)/cp)*(qi0-qi)
     qv    = qv0 + ql0 - ql + qi0 - qi

   ! Below two lines should be removed.

   ! call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
   ! qv    = max(qv,al_st*qsat+(1._r8-al_st)*1.e-12_r8) 

   ! -------------- !
   ! Send to output !
   ! -------------- !

   T_out(i)  = T
   qv_out(i) = qv
   ql_out(i) = ql
   qi_out(i) = qi
   al_st_out(i) = al_st
   ai_st_out(i) = ai_st
   ql_st_out(i) = ql_st
   qi_st_out(i) = qi_st

   enddo 

   return
   end subroutine instratus_condensate

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine instratus_core( lchnk, icol, k, p,                      &
                              T0, qv0, ql0, qi0,                      &
                              a_dc, ql_dc, qi_dc,                     & 
                              a_sc, ql_sc, qi_sc, ai_st,              &
                              qcst_crit, Tmin, Tmax, landfrac, snowh, &
                              T, qv, ql, qi )

   ! ------------------------------------------------------ !
   ! Subroutine to find saturation equilibrium state using  ! 
   ! a Newton iteration method, so that 'qc_st = qcst_crit' !
   ! is satisfied.                                          !
   ! ------------------------------------------------------ !

   use time_manager,  only: is_first_step, get_nstep

   implicit none

   integer,  intent(in)  :: lchnk      ! Chunk identifier
   integer,  intent(in)  :: icol       ! Number of atmospheric columns
   integer,  intent(in)  :: k          ! Layer index

   real(r8), intent(in)  :: p          ! Pressure [Pa]
   real(r8), intent(in)  :: T0         ! Temperature [K]
   real(r8), intent(in)  :: qv0        ! Grid-mean water vapor [kg/kg]
   real(r8), intent(in)  :: ql0        ! Grid-mean LWC [kg/kg]
   real(r8), intent(in)  :: qi0        ! Grid-mean IWC [kg/kg]

   real(r8), intent(in)  :: a_dc       ! Deep cumulus cloud fraction
   real(r8), intent(in)  :: ql_dc      ! In-deep cumulus LWC [kg/kg]
   real(r8), intent(in)  :: qi_dc      ! In-deep cumulus IWC [kg/kg]
   real(r8), intent(in)  :: a_sc       ! Shallow cumulus cloud fraction
   real(r8), intent(in)  :: ql_sc      ! In-shallow cumulus LWC [kg/kg]
   real(r8), intent(in)  :: qi_sc      ! In-shallow cumulus IWC [kg/kg]

   real(r8), intent(in)  :: ai_st      ! Ice stratus fraction (fixed)

   real(r8), intent(in)  :: Tmin       ! Minimum temperature system can have [K]
   real(r8), intent(in)  :: Tmax       ! Maximum temperature system can have [K]
   real(r8), intent(in)  :: qcst_crit  ! Critical in-stratus condensate [kg/kg]
   real(r8), intent(in)  :: landfrac   ! Land fraction
   real(r8), intent(in)  :: snowh      ! Snow depth (liquid water equivalent)

   real(r8), intent(out) :: T          ! Temperature [K]
   real(r8), intent(out) :: qv         ! Grid-mean water vapor [kg/kg]
   real(r8), intent(out) :: ql         ! Grid-mean LWC [kg/kg]
   real(r8), intent(out) :: qi         ! Grid-mean IWC [kg/kg]

   ! Local variables

   integer i                           ! Iteration index

   real(r8) muQ0, muQ
   real(r8) ql_nc0, qi_nc0, qc_nc0, qc_nc    
   real(r8) fice0, fice    
   real(r8) ficeg0, ficeg   
   real(r8) esat0, esat
   real(r8) qsat0, qsat
   real(r8) dqsdT0, dqsdT
   real(r8) dqcncdt, dastdt, dUdt
   real(r8) alpha, beta
   real(r8) U, U_nc
   real(r8) al_st_nc, G_nc
   real(r8) al_st

   ! Variables for root-finding algorithm

   integer j                          
   real(r8)  x1, x2
   real(r8)  rtsafe
   real(r8)  df, dx, dxold, f, fh, fl, temp, xh, xl
   real(r8), parameter :: xacc = 1.e-3_r8

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

   ql_nc0 = max(0._r8,ql0-a_dc*ql_dc-a_sc*ql_sc)
   qi_nc0 = max(0._r8,qi0-a_dc*qi_dc-a_sc*qi_sc)
   qc_nc0 = max(0._r8,ql0+qi0-a_dc*(ql_dc+qi_dc)-a_sc*(ql_sc+qi_sc))
   fice0  = 0._r8
   ficeg0 = 0._r8
   muQ0   = 1._r8

   ! ------------ !
   ! Root finding !
   ! ------------ !

   x1 = Tmin
   x2 = Tmax
   call funcd_instratus( x1, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
                         a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
                         qcst_crit, landfrac, snowh,                    &
                         fl, df, qc_nc, fice, al_st )
   call funcd_instratus( x2, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
                         a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
                         qcst_crit, landfrac, snowh,                    &
                         fh, df, qc_nc, fice, al_st )
   if((fl > 0._r8 .and. fh > 0._r8) .or. (fl < 0._r8 .and. fh < 0._r8)) then
       call funcd_instratus( T0, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
                             a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st, &
                             qcst_crit, landfrac, snowh,                    &
                             fl, df, qc_nc, fice, al_st )
       rtsafe = T0 
       goto 10       
   endif
   if( fl == 0._r8) then
           rtsafe = x1
           goto 10
   elseif( fh == 0._r8) then
           rtsafe = x2
           goto 10
   elseif( fl < 0._r8) then
           xl = x1
           xh = x2
   else
           xh = x1
           xl = x2
   end if
   rtsafe = 0.5_r8*(x1+x2)
   dxold = abs(x2-x1)
   dx = dxold
   call funcd_instratus( rtsafe, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
                         a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st,     &
                         qcst_crit, landfrac, snowh,                        &
                         f, df, qc_nc, fice, al_st )
   do j = 1, 20
      if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f) > 0._r8 .or. abs(2.0_r8*f) > abs(dxold*df) ) then
           dxold = dx
           dx = 0.5_r8*(xh-xl)
           rtsafe = xl + dx
           if(xl == rtsafe) goto 10
      else
           dxold = dx
           dx = f/df
           temp = rtsafe
           rtsafe = rtsafe - dx
           if (temp == rtsafe) goto 10
      end if
    ! if(abs(dx) < xacc) goto 10
      call funcd_instratus( rtsafe, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0, &
                            a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st,     &
                            qcst_crit, landfrac, snowh,                        &
                            f, df, qc_nc, fice, al_st )
    ! Sep.21.2010. Sungsu modified to enhance convergence and guarantee 'qlst_min <  qlst < qlst_max'.
      if( qcst_crit < 0.5 * ( qlst_min + qlst_max ) ) then
          if( ( qc_nc*(1._r8-fice) .gt.          qlst_min*al_st .and. &
                qc_nc*(1._r8-fice) .lt. 1.1_r8 * qlst_min*al_st ) ) goto 10
      else
          if( ( qc_nc*(1._r8-fice) .gt. 0.9_r8 * qlst_max*al_st .and. &
                qc_nc*(1._r8-fice) .lt.          qlst_max*al_st ) ) goto 10
      endif
      if(f < 0._r8) then
          xl = rtsafe
      else
          xh = rtsafe
      endif

   enddo

10 continue

   ! ------------------------------------------- !
   ! Final safety check before sending to output !
   ! ------------------------------------------- !

   qc_nc = max(0._r8,qc_nc)

   T  = rtsafe
   ql = qc_nc*(1._r8-fice) + a_dc*ql_dc + a_sc*ql_sc
   qi = qc_nc*fice + a_dc*qi_dc + a_sc*qi_sc
   qv = qv0 + ql0 + qi0 - (qc_nc + a_dc*(ql_dc+qi_dc) + a_sc*(ql_sc+qi_sc))
   qv = max(qv,1.e-12_r8) 

   return
   end subroutine instratus_core

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine funcd_instratus( T, p, T0, qv0, ql0, qi0, fice0, muQ0, qc_nc0,   &
                               a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, ai_st,  &
                               qcst_crit, landfrac, snowh,                     &
                               f, fg, qc_nc, fice, al_st ) 

   ! --------------------------------------------------- !
   ! Subroutine to find function value and gradient at T !
   ! --------------------------------------------------- !

   implicit none

   real(r8), intent(in)  :: T          ! Iteration temperature [K]

   real(r8), intent(in)  :: p          ! Pressure [Pa]
   real(r8), intent(in)  :: T0         ! Initial temperature [K]
   real(r8), intent(in)  :: qv0        ! Grid-mean water vapor [kg/kg]
   real(r8), intent(in)  :: ql0        ! Grid-mean LWC [kg/kg]
   real(r8), intent(in)  :: qi0        ! Grid-mean IWC [kg/kg]
   real(r8), intent(in)  :: fice0      ! 
   real(r8), intent(in)  :: muQ0       ! 
   real(r8), intent(in)  :: qc_nc0     ! 

   real(r8), intent(in)  :: a_dc       ! Deep cumulus cloud fraction
   real(r8), intent(in)  :: ql_dc      ! In-deep cumulus LWC [kg/kg]
   real(r8), intent(in)  :: qi_dc      ! In-deep cumulus IWC [kg/kg]
   real(r8), intent(in)  :: a_sc       ! Shallow cumulus cloud fraction
   real(r8), intent(in)  :: ql_sc      ! In-shallow cumulus LWC [kg/kg]
   real(r8), intent(in)  :: qi_sc      ! In-shallow cumulus IWC [kg/kg]

   real(r8), intent(in)  :: ai_st      ! Ice stratus fraction (fixed)

   real(r8), intent(in)  :: qcst_crit  ! Critical in-stratus condensate [kg/kg]
   real(r8), intent(in)  :: landfrac   ! Land fraction
   real(r8), intent(in)  :: snowh      ! Snow depth (liquid water equivalent)

   real(r8), intent(out) :: f          ! Value of minimization function at T
   real(r8), intent(out) :: fg         ! Gradient of minimization function 
   real(r8), intent(out) :: qc_nc      !
   real(r8), intent(out) :: al_st      !
   real(r8), intent(out) :: fice       !

   ! Local variables

   real(r8) esat
   real(r8) qsat
   real(r8) dqsdT
   real(r8) dqcncdt
   real(r8) alpha
   real(r8) beta
   real(r8) U
   real(r8) U_nc
   real(r8) al_st_nc
   real(r8) G_nc
   real(r8) dUdt
   real(r8) dalstdt
   real(r8) qv

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

   call vqsatd2_water_single(T,p,esat,qsat,dqsdT)

   fice    = fice0 
   qc_nc   = (cp/hlatv)*(T-T0)+muQ0*qc_nc0       
   dqcncdt = (cp/hlatv) 
   qv      = (qv0 + ql0 + qi0 - (qc_nc + a_dc*(ql_dc+qi_dc) + a_sc*(ql_sc+qi_sc)))
   alpha   = (1._r8/qsat)
   beta    = (qv/qsat**2._r8)*dqsdT 

   U      =  (qv/qsat)
   U_nc   =   U
   if( CAMstfrac ) then
       call astG_RHU_single(U_nc,p,qv,landfrac,snowh,al_st_nc,G_nc)
   else
       call astG_PDF_single(U_nc,p,qv,landfrac,snowh,al_st_nc,G_nc)
   endif
   al_st   =  (1._r8-a_dc-a_sc)*al_st_nc 
   dUdt    = -(alpha*dqcncdt+beta)
   dalstdt =  (1._r8/G_nc)*dUdt
   if( U_nc .eq. 1._r8 ) dalstdt = 0._r8

   f  = qc_nc   - qcst_crit*al_st
   fg = dqcncdt - qcst_crit*dalstdt

   return
   end subroutine funcd_instratus

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine gridmean_RH( lchnk, icol, k, p, T, qv, ql, qi,       &
                           a_dc, ql_dc, qi_dc, a_sc, ql_sc, qi_sc, &
                           landfrac, snowh )

   ! ------------------------------------------------------------- !
   ! Subroutine to force grid-mean RH = 1 when RH > 1              !
   ! This is condensation process similar to instratus_condensate. !
   ! During condensation, we assume 'fice' is maintained in this   !
   ! verison for MG not for RK.                                    !
   ! ------------------------------------------------------------- !

   use time_manager,  only: is_first_step, get_nstep

   implicit none

   integer,  intent(in)    :: lchnk      ! Chunk identifier
   integer,  intent(in)    :: icol       ! Number of atmospheric columns
   integer,  intent(in)    :: k          ! Layer index

   real(r8), intent(in)    :: p          ! Pressure [Pa]
   real(r8), intent(inout) :: T          ! Temperature [K]
   real(r8), intent(inout) :: qv         ! Grid-mean water vapor [kg/kg]
   real(r8), intent(inout) :: ql         ! Grid-mean LWC [kg/kg]
   real(r8), intent(inout) :: qi         ! Grid-mean IWC [kg/kg]

   real(r8), intent(in)    :: a_dc       ! Deep cumulus cloud fraction
   real(r8), intent(in)    :: ql_dc      ! In-deep cumulus LWC [kg/kg]
   real(r8), intent(in)    :: qi_dc      ! In-deep cumulus IWC [kg/kg]
   real(r8), intent(in)    :: a_sc       ! Shallow cumulus cloud fraction
   real(r8), intent(in)    :: ql_sc      ! In-shallow cumulus LWC [kg/kg]
   real(r8), intent(in)    :: qi_sc      ! In-shallow cumulus IWC [kg/kg]

   real(r8), intent(in)    :: landfrac   ! Land fraction
   real(r8), intent(in)    :: snowh      ! Snow depth (liquid water equivalent)

   ! Local variables

   integer m                             ! Iteration index

   real(r8)  ql_nc0, qi_nc0, qc_nc0
   real(r8)  Tscale
   real(r8)  Tc, qt, qc, dqcdt, qc_nc    
   real(r8)  esat, qsat, dqsdT
   real(r8)  al_st_nc, G_nc
   real(r8)  f, fg
   real(r8), parameter :: xacc = 1.e-3_r8

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

   ql_nc0 = max(0._r8,ql-a_dc*ql_dc-a_sc*ql_sc)
   qi_nc0 = max(0._r8,qi-a_dc*qi_dc-a_sc*qi_sc)
   qc_nc0 = max(0._r8,ql+qi-a_dc*(ql_dc+qi_dc)-a_sc*(ql_sc+qi_sc))
   Tc    = T - (hlatv/cp)*ql
   qt    = qv + ql

   do m = 1, 20
      call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
      Tscale = hlatv/cp
      qc     = (T-Tc)/Tscale
      dqcdt  = 1._r8/Tscale
      f      = qsat + qc - qt 
      fg     = dqsdT + dqcdt
      fg     = sign(1._r8,fg)*max(1.e-10_r8,abs(fg))
    ! Sungsu modified convergence criteria to speed up convergence and guarantee RH <= 1.
      if( qc .ge. 0._r8 .and. ( qt - qc ) .ge. 0.999_r8*qsat .and. ( qt - qc ) .le. 1._r8*qsat ) then
          goto 10
      endif
      T = T - f/fg
   enddo
 ! write(iulog,*) 'Convergence in gridmean_RH is not reached. RH = ', ( qt - qc ) / qsat
10 continue

   call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
 ! Sungsu modified 'qv = qsat' in consistent with the modified convergence criteria above.
   qv = min(qt,qsat) ! Modified
   ql = qt - qv
   T  = Tc + (hlatv/cp)*ql

   return
   end subroutine gridmean_RH

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine positive_moisture( ncol, dt, qvmin, qlmin, qimin, dp, qv, ql, qi, t, qvten, qlten, qiten, tten )

   ! ------------------------------------------------------------------------------- !
   ! If any 'ql < qlmin, qi < qimin, qv < qvmin' are developed in any layer,         !
   ! force them to be larger than minimum value by (1) condensating water vapor      !
   ! into liquid or ice, and (2) by transporting water vapor from the very lower     !
   ! layer. '2._r8' is multiplied to the minimum values for safety.                  !
   ! Update final state variables and tendencies associated with this correction.    !
   ! If any condensation happens, update (s,t) too.                                  !
   ! Note that (qv,ql,qi,t,s) are final state variables after applying corresponding !
   ! input tendencies.                                                               !
   ! Be careful the order of k : '1': top layer, 'pver' : near-surface layer         ! 
   ! ------------------------------------------------------------------------------- !

   implicit none
   integer,  intent(in)     :: ncol
   real(r8), intent(in)     :: dt
   real(r8), intent(in)     :: dp(pcols,pver), qvmin(pcols,pver), qlmin(pcols,pver), qimin(pcols,pver)
   real(r8), intent(inout)  :: qv(pcols,pver), ql(pcols,pver), qi(pcols,pver), t(pcols,pver)
   real(r8), intent(out)    :: qvten(pcols,pver), qlten(pcols,pver), qiten(pcols,pver), tten(pcols,pver)
   integer   i, k
   real(r8)  dql, dqi, dqv, sum, aa, dum 

   tten(:ncol,:pver)  = 0._r8
   qvten(:ncol,:pver) = 0._r8
   qlten(:ncol,:pver) = 0._r8
   qiten(:ncol,:pver) = 0._r8

   do i = 1, ncol
      do k = 1, pver
         if( qv(i,k) .lt. qvmin(i,k) .or. ql(i,k) .lt. qlmin(i,k) .or. qi(i,k) .lt. qimin(i,k) ) then
             goto 10
         endif
      enddo
      goto 11
   10 continue
      do k = 1, pver    ! From the top to the 1st (lowest) layer from the surface
         dql = max(0._r8,1._r8*qlmin(i,k)-ql(i,k))
         dqi = max(0._r8,1._r8*qimin(i,k)-qi(i,k))
         qlten(i,k) = qlten(i,k) +  dql/dt
         qiten(i,k) = qiten(i,k) +  dqi/dt
         qvten(i,k) = qvten(i,k) - (dql+dqi)/dt
         tten(i,k)  = tten(i,k)  + (hlatv/cp)*(dql/dt) + ((hlatv+hlatf)/cp)*(dqi/dt)
         ql(i,k)    = ql(i,k) + dql
         qi(i,k)    = qi(i,k) + dqi
         qv(i,k)    = qv(i,k) - dql - dqi
         t(i,k)     = t(i,k)  + (hlatv * dql + (hlatv+hlatf) * dqi)/cp
         dqv        = max(0._r8,1._r8*qvmin(i,k)-qv(i,k))
         qvten(i,k) = qvten(i,k) + dqv/dt
         qv(i,k)    = qv(i,k)    + dqv
         if( k .ne. pver ) then 
             qv(i,k+1)    = qv(i,k+1)    - dqv*dp(i,k)/dp(i,k+1)
             qvten(i,k+1) = qvten(i,k+1) - dqv*dp(i,k)/dp(i,k+1)/dt
         endif
         qv(i,k) = max(qv(i,k),qvmin(i,k))
         ql(i,k) = max(ql(i,k),qlmin(i,k))
         qi(i,k) = max(qi(i,k),qimin(i,k))
      end do
      ! Extra moisture used to satisfy 'qv(i,pver)=qvmin' is proportionally 
      ! extracted from all the layers that has 'qv > 2*qvmin'. This fully
      ! preserves column moisture. 
      if( dqv .gt. 1.e-20 ) then
          sum = 0._r8
          do k = 1, pver
             if( qv(i,k) .gt. 2._r8*qvmin(i,k) ) sum = sum + qv(i,k)*dp(i,k)
          enddo
          aa = dqv*dp(i,pver)/max(1.e-20_r8,sum)
          if( aa .lt. 0.5_r8 ) then
              do k = 1, pver
                 if( qv(i,k) .gt. 2._r8*qvmin(i,k) ) then
                     dum        = aa*qv(i,k)
                     qv(i,k)    = qv(i,k) - dum
                     qvten(i,k) = qvten(i,k) - dum/dt
                 endif
              enddo 
          else 
              write(iulog,*) 'Full positive_moisture is impossible in Park Macro'
          endif
      endif 
11 continue
   enddo
   return

   end subroutine positive_moisture

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine astG_PDF_single( U, p, qv, landfrac, snowh, a, Ga )

   ! --------------------------------------------------------- !
   ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the     !
   ! analytical formulation of triangular PDF.                 !
   ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', !
   ! so using constant 'dV' assume that width is proportional  !
   ! to the saturation specific humidity.                      !
   !    dV ~ 0.1.                                              !
   !    cldrh : RH of in-stratus( = 1 if no supersaturation)   !
   ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
   ! G is discontinuous across U = 1.  In fact, it does not    !
   ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that   !
   ! they will produce the same results.                       !
   ! --------------------------------------------------------- !

   implicit none

   real(r8), intent(in)  :: U                     ! Relative humidity
   real(r8), intent(in)  :: p                     ! Pressure [Pa]
   real(r8), intent(in)  :: qv                    ! Grid-mean water vapor specific humidity [kg/kg]
   real(r8), intent(in)  :: landfrac              ! Land fraction
   real(r8), intent(in)  :: snowh                 ! Snow depth (liquid water equivalent)

   real(r8), intent(out) :: a                     ! Stratus fraction
   real(r8), intent(out) :: Ga                    ! dU/da

   ! Local variables
   integer :: i                                   ! Loop indexes
   real(r8) dV                                    ! Width of triangular PDF
   real(r8) cldrh                                 ! RH of stratus cloud
   real(r8) rhmin                                 ! Critical RH
   real(r8) rhwght
                            
   ! Statement functions
   logical land
   land = nint(landfrac) == 1

   ! ---------- !
   ! Parameters !
   ! ---------- !

   cldrh  = 1.0_r8

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

   if( p .ge. premib ) then

       if( land .and. (snowh.le.0.000001_r8) ) then
           rhmin = rhminl - 0.10_r8
       else
           rhmin = rhminl
       endif

       dV = cldrh - rhmin

       if( U .ge. 1._r8 ) then
           a  = 1._r8
           Ga = 1.e10_r8
       elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
           a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
           Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
       elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
           a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & 
                      (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
           Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
       elseif( U .le. (cldrh-dV) ) then
           a  = 0._r8
           Ga = 1.e10_r8
       endif

       if( freeze_dry ) then
           a  = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
           Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) 
       endif

   elseif( p .lt. premit ) then

       rhmin = rhminh
       dV    = cldrh - rhmin

       if( U .ge. 1._r8 ) then
           a  = 1._r8
           Ga = 1.e10_r8
       elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
           a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
           Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
       elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
           a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & 
                      (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
           Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
       elseif( U .le. (cldrh-dV) ) then
           a  = 0._r8
           Ga = 1.e10_r8
       endif

   else

       rhwght = (premib-(max(p,premit)))/(premib-premit)

     ! if( land .and. (snowh.le.0.000001_r8) ) then
     !     rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
     ! else
           rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
     ! endif

       dV    = cldrh - rhmin

       if( U .ge. 1._r8 ) then
           a  = 1._r8
           Ga = 1.e10_r8
       elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
           a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
           Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
       elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
           a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & 
                         (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
           Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
       elseif( U .le. (cldrh-dV) ) then
           a  = 0._r8
           Ga = 1.e10_r8
       endif

   endif

   return
   end subroutine astG_PDF_single

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine astG_PDF( U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol )

   ! --------------------------------------------------------- !
   ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the     !
   ! analytical formulation of triangular PDF.                 !
   ! Here, 'dV' is the ratio of 'half-width of PDF / qs(p,T)', !
   ! so using constant 'dV' assume that width is proportional  !
   ! to the saturation specific humidity.                      !
   !    dV ~ 0.1.                                              !
   !    cldrh : RH of in-stratus( = 1 if no supersaturation)   !
   ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
   ! G is discontinuous across U = 1.  In fact, it does not    !
   ! matter whether Ga = 1.e10 or 0 at a = 1: I derived that   !
   ! they will produce the same results.                       !
   ! --------------------------------------------------------- !

   implicit none

   integer,  intent(in)  :: ncol
   real(r8), intent(in)  :: U_in(pcols)           ! Relative humidity
   real(r8), intent(in)  :: p_in(pcols)           ! Pressure [Pa]
   real(r8), intent(in)  :: qv_in(pcols)          ! Grid-mean water vapor specific humidity [kg/kg]
   real(r8), intent(in)  :: landfrac_in(pcols)    ! Land fraction
   real(r8), intent(in)  :: snowh_in(pcols)       ! Snow depth (liquid water equivalent)

   real(r8), intent(out) :: a_out(pcols)          ! Stratus fraction
   real(r8), intent(out) :: Ga_out(pcols)         ! dU/da

   real(r8)              :: U                     ! Relative humidity
   real(r8)              :: p                     ! Pressure [Pa]
   real(r8)              :: qv                    ! Grid-mean water vapor specific humidity [kg/kg]
   real(r8)              :: landfrac              ! Land fraction
   real(r8)              :: snowh                 ! Snow depth (liquid water equivalent)

   real(r8)              :: a                     ! Stratus fraction
   real(r8)              :: Ga                    ! dU/da

   ! Local variables
   integer :: i                                   ! Loop indexes
   real(r8) dV                                    ! Width of triangular PDF
   real(r8) cldrh                                 ! RH of stratus cloud
   real(r8) rhmin                                 ! Critical RH
   real(r8) rhwght
                            
   ! Statement functions
   logical land
   land(i) = nint(landfrac_in(i)) == 1

   ! ---------- !
   ! Parameters !
   ! ---------- !

   cldrh  = 1.0_r8

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

   a_out(:)  = 0._r8
   Ga_out(:) = 0._r8

   do i = 1, ncol

   U        = U_in(i)      
   p        = p_in(i)        
   qv       = qv_in(i)       
   landfrac = landfrac_in(i) 
   snowh    = snowh_in(i)    

   if( p .ge. premib ) then

       if( land(i) .and. (snowh.le.0.000001_r8) ) then
           rhmin = rhminl - 0.10_r8
       else
           rhmin = rhminl
       endif

       dV = cldrh - rhmin

       if( U .ge. 1._r8 ) then
           a  = 1._r8
           Ga = 1.e10_r8
       elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
           a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
           Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
       elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
           a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & 
                      (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
           Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
       elseif( U .le. (cldrh-dV) ) then
           a  = 0._r8
           Ga = 1.e10_r8
       endif

       if( freeze_dry ) then
           a  = a *max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
           Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) 
       endif

   elseif( p .lt. premit ) then

       rhmin = rhminh
       dV    = cldrh - rhmin

       if( U .ge. 1._r8 ) then
           a  = 1._r8
           Ga = 1.e10_r8
       elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
           a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
           Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
       elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
           a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & 
                      (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
           Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
       elseif( U .le. (cldrh-dV) ) then
           a  = 0._r8
           Ga = 1.e10_r8
       endif

   else

       rhwght = (premib-(max(p,premit)))/(premib-premit)

     ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
     !     rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
     ! else
           rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
     ! endif

       dV    = cldrh - rhmin

       if( U .ge. 1._r8 ) then
           a  = 1._r8
           Ga = 1.e10_r8
       elseif( U .gt. (cldrh-dV/6._r8) .and. U .lt. 1._r8 ) then
           a  = 1._r8 - (-3._r8/sqrt(2._r8)*(U-cldrh)/dV)**(2._r8/3._r8)
           Ga = dV/sqrt(2._r8)*sqrt(1._r8-a)
       elseif( U .gt. (cldrh-dV) .and. U .le. (cldrh-dV/6._r8) ) then
           a  = 4._r8*(cos((1._r8/3._r8)*(acos((3._r8/2._r8/sqrt(2._r8))* & 
                         (1._r8+(U-cldrh)/dV))-2._r8*3.141592_r8)))**2._r8
           Ga = dV/sqrt(2._r8)*(1._r8/sqrt(a)-sqrt(a))
       elseif( U .le. (cldrh-dV) ) then
           a  = 0._r8
           Ga = 1.e10_r8
       endif

   endif

   a_out(i)  = a
   Ga_out(i) = Ga 

   enddo

   return
   end subroutine astG_PDF

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine astG_RHU_single( U, p, qv, landfrac, snowh, a, Ga )

   ! --------------------------------------------------------- !
   ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the     !
   ! CAM35 cloud fraction formula.                             !
   ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core !  
   ! For the other cases, I should re-define 'rhminl,rhminh' & !
   ! 'premib,premit'.                                          !
   ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
   ! G is discontinuous across U = 1.                          !
   ! --------------------------------------------------------- !

   implicit none

   real(r8), intent(in)  :: U               ! Relative humidity
   real(r8), intent(in)  :: p               ! Pressure [Pa]
   real(r8), intent(in)  :: qv              ! Grid-mean water vapor specific humidity [kg/kg]
   real(r8), intent(in)  :: landfrac        ! Land fraction
   real(r8), intent(in)  :: snowh           ! Snow depth (liquid water equivalent)

   real(r8), intent(out) :: a               ! Stratus fraction
   real(r8), intent(out) :: Ga              ! dU/da

   ! Local variables
   real(r8) rhmin                                 ! Critical RH
   real(r8) rhdif                                 ! Factor for stratus fraction
   real(r8) rhwght

   ! Statement functions
   logical land
   land = nint(landfrac) == 1

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

   if( p .ge. premib ) then

       if( land .and. (snowh.le.0.000001_r8) ) then
           rhmin = rhminl - 0.10_r8
       else
           rhmin = rhminl
       endif
       rhdif = (U-rhmin)/(1.0_r8-rhmin)
       a  = min(1._r8,(max(rhdif,0.0_r8))**2) 
       if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
            Ga = 1.e20_r8
       else          
            Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
       endif
       if( freeze_dry ) then
           a  = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
           Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) 
       endif

   elseif( p .lt. premit ) then

       rhmin = rhminh
       rhdif = (U-rhmin)/(1.0_r8-rhmin)
       a  = min(1._r8,(max(rhdif,0._r8))**2)
       if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
            Ga = 1.e20_r8
       else          
            Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
       endif

   else

       rhwght = (premib-(max(p,premit)))/(premib-premit)

     ! if( land .and. (snowh.le.0.000001_r8) ) then
     !     rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
     ! else
           rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
     ! endif

       rhdif = (U-rhmin)/(1.0_r8-rhmin)
       a  = min(1._r8,(max(rhdif,0._r8))**2)
       if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
            Ga = 1.e10_r8
       else          
            Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
       endif

   endif

   return
   end subroutine astG_RHU_single

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine astG_RHU( U_in, p_in, qv_in, landfrac_in, snowh_in, a_out, Ga_out, ncol )

   ! --------------------------------------------------------- !
   ! Compute 'stratus fraction(a)' and Gs=(dU/da) from the     !
   ! CAM35 cloud fraction formula.                             !
   ! Below is valid only for CAMUW at 1.9x2.5 fv dynamics core !  
   ! For the other cases, I should re-define 'rhminl,rhminh' & !
   ! 'premib,premit'.                                          !
   ! Note that if U > 1, Ga = 1.e10 instead of Ga = 0, that is !
   ! G is discontinuous across U = 1.                          !
   ! --------------------------------------------------------- !

   implicit none

   integer,  intent(in)  :: ncol
   real(r8), intent(in)  :: U_in(pcols)           ! Relative humidity
   real(r8), intent(in)  :: p_in(pcols)           ! Pressure [Pa]
   real(r8), intent(in)  :: qv_in(pcols)          ! Grid-mean water vapor specific humidity [kg/kg]
   real(r8), intent(in)  :: landfrac_in(pcols)    ! Land fraction
   real(r8), intent(in)  :: snowh_in(pcols)       ! Snow depth (liquid water equivalent)

   real(r8), intent(out) :: a_out(pcols)          ! Stratus fraction
   real(r8), intent(out) :: Ga_out(pcols)         ! dU/da

   real(r8)              :: U                     ! Relative humidity
   real(r8)              :: p                     ! Pressure [Pa]
   real(r8)              :: qv                    ! Grid-mean water vapor specific humidity [kg/kg]
   real(r8)              :: landfrac              ! Land fraction
   real(r8)              :: snowh                 ! Snow depth (liquid water equivalent)

   real(r8)              :: a                     ! Stratus fraction
   real(r8)              :: Ga                    ! dU/da

   ! Local variables
   integer  i
   real(r8) rhmin                                 ! Critical RH
   real(r8) rhdif                                 ! Factor for stratus fraction
   real(r8) rhwght

   ! Statement functions
   logical land
   land(i) = nint(landfrac_in(i)) == 1

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

   a_out(:) = 0._r8
   Ga_out(:) = 0._r8

   do i = 1, ncol

   U        = U_in(i)      
   p        = p_in(i)        
   qv       = qv_in(i)       
   landfrac = landfrac_in(i) 
   snowh    = snowh_in(i)    

   if( p .ge. premib ) then

       if( land(i) .and. (snowh.le.0.000001_r8) ) then
           rhmin = rhminl - 0.10_r8
       else
           rhmin = rhminl
       endif
       rhdif = (U-rhmin)/(1.0_r8-rhmin)
       a  = min(1._r8,(max(rhdif,0.0_r8))**2) 
       if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
            Ga = 1.e20_r8
       else          
            Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
       endif
       if( freeze_dry ) then
           a  = a*max(0.15_r8,min(1.0_r8,qv/0.0030_r8))
           Ga = Ga/max(0.15_r8,min(1.0_r8,qv/0.0030_r8)) 
       endif

   elseif( p .lt. premit ) then

       rhmin = rhminh
       rhdif = (U-rhmin)/(1.0_r8-rhmin)
       a  = min(1._r8,(max(rhdif,0._r8))**2)
       if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
            Ga = 1.e20_r8
       else          
            Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
       endif

   else

       rhwght = (premib-(max(p,premit)))/(premib-premit)

     ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
     !     rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
     ! else
           rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
     ! endif

       rhdif = (U-rhmin)/(1.0_r8-rhmin)
       a  = min(1._r8,(max(rhdif,0._r8))**2)
       if( (U.ge.1._r8) .or. (U.le.rhmin) ) then
            Ga = 1.e10_r8
       else          
            Ga = 0.5_r8*(1._r8-rhmin)*((1._r8-rhmin)/(U-rhmin))
       endif

   endif

   a_out(i)  = a
   Ga_out(i) = Ga 

   enddo

   return
   end subroutine astG_RHU

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine aist_single( qv, T, p, qi, landfrac, snowh, aist )

   ! --------------------------------------------------------- !
   ! Compute non-physical ice stratus fraction                 ! 
   ! --------------------------------------------------------- !

   use physconst,     only: rair
   use wv_saturation, only: vqsatd2_water_single

   implicit none
  
   real(r8), intent(in)  :: qv              ! Grid-mean water vapor[kg/kg]
   real(r8), intent(in)  :: T               ! Temperature
   real(r8), intent(in)  :: p               ! Pressure [Pa]
   real(r8), intent(in)  :: qi              ! Grid-mean ice water content [kg/kg]
   real(r8), intent(in)  :: landfrac        ! Land fraction
   real(r8), intent(in)  :: snowh           ! Snow depth (liquid water equivalent)

   real(r8), intent(out) :: aist            ! Non-physical ice stratus fraction ( 0<= aist <= 1 )

   ! Local variables
   real(r8) rhmin                           ! Critical RH
   real(r8) rhwght

   real(r8) a,b,c,as,bs,cs                  ! Fit parameters
   real(r8) Kc                              ! Constant for ice cloud calc (wood & field)
   real(r8) ttmp                            ! Limited temperature
   real(r8) icicval                         ! Empirical IWC value [ kg/kg ]
   real(r8) rho                             ! Local air density
   real(r8) esl                             ! Liq sat vapor pressure
   real(r8) esi                             ! Ice sat vapor pressure
   real(r8) ncf,phi                         ! Wilson and Ballard parameters
   real(r8) esat, qsat, dqsdT

   real(r8) rhi                             ! grid box averaged relative humidity over ice
   real(r8) minice                          ! minimum grid box avg ice for having a 'cloud'
   real(r8) mincld                          ! minimum ice cloud fraction threshold
   real(r8) icimr                           ! in cloud ice mixing ratio
 ! real(r8) qist_min                        ! minimum in cloud ice mixing ratio
 ! real(r8) qist_max                        ! maximum in cloud ice mixing ratio                
   real(r8) rhdif                           ! working variable for slingo scheme

 ! integer iceopt                           ! option for ice cloud closure 
                                            ! 1=wang & sassen 2=schiller (iciwc)  
                                            ! 3=wood & field, 4=Wilson (based on smith)
                                            ! 5=modified slingo (ssat & empyt cloud)

   real(r8) icecrit                         ! Critical RH for ice clouds in Wilson & Ballard closure ( smaller = more ice clouds )

   ! Statement functions
   logical land
   land = nint(landfrac) == 1

   ! --------- !
   ! Constants !
   ! --------- !

   ! iceopt = 5

   ! Wang and Sassen IWC paramters ( Option.1 )
     a = 26.87_r8
     b = 0.569_r8
     c = 0.002892_r8
   ! Schiller parameters ( Option.2 )
     as = -68.4202_r8
     bs = 0.983917_r8
     cs = 2.81795_r8
   ! Wood and Field parameters ( Option.3 )
     Kc = 75._r8
   ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds)
     icecrit = 0.93_r8
   ! Slingo modified (option 5)
     minice = 1.e-12_r8
     mincld = 1.e-4_r8
   ! qist_min = 1.e-7_r8 
   ! qist_max = 5.e-3_r8

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

     call vqsatd2_water_single(T,p,esat,qsat,dqsdT)
     esl = polysvp(T,0)
     esi = polysvp(T,1)
          
     if( iceopt.lt.3 ) then
         if( iceopt.eq.1 ) then
             ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8
             icicval = a + b * ttmp + c * ttmp**2._r8
             rho = p/(rair*T)
             icicval = icicval * 1.e-6_r8 / rho 
         else
             ttmp = max(190._r8,min(T,273.16_r8))
             icicval = 10._r8 **(as * bs**ttmp + cs)
             icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
         endif
         aist =  max(0._r8,min(qi/icicval,1._r8)) 
     elseif( iceopt.eq.3 ) then
         aist = 1._r8 - exp(-Kc*qi/(qsat*(esi/esl)))
         aist = max(0._r8,min(aist,1._r8))
     elseif( iceopt.eq.4) then
         if( p .ge. premib ) then
             if( land .and. (snowh.le.0.000001_r8) ) then
                 rhmin = rhminl - 0.10_r8
             else
                 rhmin = rhminl
             endif
         elseif( p .lt. premit ) then
             rhmin = rhminh
         else
             rhwght = (premib-(max(p,premit)))/(premib-premit)
           ! if( land .and. (snowh.le.0.000001_r8) ) then
           !     rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
           ! else
                 rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
           ! endif
         endif
         ncf = qi/((1._r8 - icecrit)*qsat)
         if( ncf.le.0._r8 ) then 
             aist = 0._r8
         elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then 
             aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
         elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then
             phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
             aist = (1._r8 - 4._r8 * cos(phi) * cos(phi))
         else
             aist = 1._r8
         endif
             aist = max(0._r8,min(aist,1._r8))
     elseif (iceopt.eq.5) then 
! set rh ice cloud fraction
             rhi= (qv+qi)/qsat * (esl/esi)
             rhdif= (rhi-rhmini) / (rhmaxi - rhmini)
             aist = min(1.0_r8, max(rhdif,0._r8)**2)

! limiter to remove empty cloud and ice with no cloud
! and set icecld fraction to mincld if ice exists

             if (qi.lt.minice) then
                aist=0._r8
             else
                aist=max(mincld,aist)
             endif

! enforce limits on icimr
             if (qi.ge.minice) then
                icimr=qi/aist

!minimum
                if (icimr.lt.qist_min) then
                   aist = max(0._r8,min(1._r8,qi/qist_min))
                endif
!maximum
                if (icimr.gt.qist_max) then
                   aist = max(0._r8,min(1._r8,qi/qist_max))
                endif

             endif
     endif 

   ! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate
   ! computed after updating 'qi_st'.  

     aist = max(0._r8,min(aist,0.999_r8))

   return
   end subroutine aist_single

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

   subroutine aist_vector( qv_in, T_in, p_in, qi_in, landfrac_in, snowh_in, aist_out, ncol )

   ! --------------------------------------------------------- !
   ! Compute non-physical ice stratus fraction                 ! 
   ! --------------------------------------------------------- !

   use physconst,     only: rair
   use wv_saturation, only: vqsatd2_water

   implicit none
  
   integer,  intent(in)  :: ncol 
   real(r8), intent(in)  :: qv_in(pcols)       ! Grid-mean water vapor[kg/kg]
   real(r8), intent(in)  :: T_in(pcols)        ! Temperature
   real(r8), intent(in)  :: p_in(pcols)        ! Pressure [Pa]
   real(r8), intent(in)  :: qi_in(pcols)       ! Grid-mean ice water content [kg/kg]
   real(r8), intent(in)  :: landfrac_in(pcols) ! Land fraction
   real(r8), intent(in)  :: snowh_in(pcols)    ! Snow depth (liquid water equivalent)
   real(r8), intent(out) :: aist_out(pcols)    ! Non-physical ice stratus fraction ( 0<= aist <= 1 )

   ! Local variables

   real(r8) qv                              ! Grid-mean water vapor[kg/kg]
   real(r8) T                               ! Temperature
   real(r8) p                               ! Pressure [Pa]
   real(r8) qi                              ! Grid-mean ice water content [kg/kg]
   real(r8) landfrac                        ! Land fraction
   real(r8) snowh                           ! Snow depth (liquid water equivalent)
   real(r8) aist                            ! Non-physical ice stratus fraction ( 0<= aist <= 1 )

   real(r8) rhmin                           ! Critical RH
   real(r8) rhwght

   real(r8) a,b,c,as,bs,cs                  ! Fit parameters
   real(r8) Kc                              ! Constant for ice cloud calc (wood & field)
   real(r8) ttmp                            ! Limited temperature
   real(r8) icicval                         ! Empirical IWC value [ kg/kg ]
   real(r8) rho                             ! Local air density
   real(r8) esl                             ! Liq sat vapor pressure
   real(r8) esi                             ! Ice sat vapor pressure
   real(r8) ncf,phi                         ! Wilson and Ballard parameters
   real(r8) esat, qsat, dqsdT
   real(r8) esat_in(pcols)
   real(r8) qsat_in(pcols)
   real(r8) dqsdT_in(pcols)

   real(r8) rhi                             ! grid box averaged relative humidity over ice
   real(r8) minice                          ! minimum grid box avg ice for having a 'cloud'
   real(r8) mincld                          ! minimum ice cloud fraction threshold
   real(r8) icimr                           ! in cloud ice mixing ratio
 ! real(r8) qist_min                        ! minimum in cloud ice mixing ratio
 ! real(r8) qist_max                        ! maximum in cloud ice mixing ratio                
   real(r8) rhdif                           ! working variable for slingo scheme

   integer i
 ! integer iceopt                           ! option for ice cloud closure 
                                            ! 1=wang & sassen 2=schiller (iciwc)  
                                            ! 3=wood & field, 4=Wilson (based on smith)
                                            ! 5=modified slingo (ssat & empyt cloud)

   real(r8) icecrit                         ! Critical RH for ice clouds in Wilson & Ballard closure ( smaller = more ice clouds )

   ! Statement functions
   logical land
   land(i) = nint(landfrac_in(i)) == 1

   ! --------- !
   ! Constants !
   ! --------- !

   ! iceopt = 5

   ! Wang and Sassen IWC paramters ( Option.1 )
     a = 26.87_r8
     b = 0.569_r8
     c = 0.002892_r8
   ! Schiller parameters ( Option.2 )
     as = -68.4202_r8
     bs = 0.983917_r8
     cs = 2.81795_r8
   ! Wood and Field parameters ( Option.3 )
     Kc = 75._r8
   ! Wilson & Ballard closure ( Option.4. smaller = more ice clouds)
     icecrit = 0.93_r8
   ! Slingo modified (option 5)
     minice = 1.e-12_r8
     mincld = 1.e-4_r8
   ! qist_min = 1.e-7_r8
   ! qist_max = 5.e-3_r8

   ! ---------------- !
   ! Main computation !
   ! ---------------- !

     aist_out(:) = 0._r8
     esat_in(:)  = 0._r8
     qsat_in(:)  = 0._r8
     dqsdT_in(:) = 0._r8

     call vqsatd2_water(T_in(1:ncol),p_in(1:ncol),esat_in(1:ncol),qsat_in(1:ncol),dqsdT_in(1:ncol),ncol)
     
     do i = 1, ncol

     landfrac = landfrac_in(i)     
     snowh = snowh_in(i)   
     T = T_in(i)
     qv = qv_in(i)
     p = p_in(i)
     qi = qi_in(i)
     qsat = qsat_in(i)
     esl = polysvp(T,0)
     esi = polysvp(T,1)
          
     if( iceopt.lt.3 ) then
         if( iceopt.eq.1 ) then
             ttmp = max(195._r8,min(T,253._r8)) - 273.16_r8
             icicval = a + b * ttmp + c * ttmp**2._r8
             rho = p/(rair*T)
             icicval = icicval * 1.e-6_r8 / rho 
         else
             ttmp = max(190._r8,min(T,273.16_r8))
             icicval = 10._r8 **(as * bs**ttmp + cs)
             icicval = icicval * 1.e-6_r8 * 18._r8 / 28.97_r8
         endif
         aist =  max(0._r8,min(qi/icicval,1._r8)) 
     elseif( iceopt.eq.3 ) then
         aist = 1._r8 - exp(-Kc*qi/(qsat*(esi/esl)))
         aist = max(0._r8,min(aist,1._r8))
     elseif( iceopt.eq.4) then
         if( p .ge. premib ) then
             if( land(i) .and. (snowh.le.0.000001_r8) ) then
                 rhmin = rhminl - 0.10_r8
             else
                 rhmin = rhminl
             endif
         elseif( p .lt. premit ) then
             rhmin = rhminh
         else
             rhwght = (premib-(max(p,premit)))/(premib-premit)
           ! if( land(i) .and. (snowh.le.0.000001_r8) ) then
           !     rhmin = rhminh*rhwght + (rhminl - 0.10_r8)*(1.0_r8-rhwght)
           ! else
                 rhmin = rhminh*rhwght + rhminl*(1.0_r8-rhwght)
           ! endif
         endif
         ncf = qi/((1._r8 - icecrit)*qsat)
         if( ncf.le.0._r8 ) then 
             aist = 0._r8
         elseif( ncf.gt.0._r8 .and. ncf.le.1._r8/6._r8 ) then 
             aist = 0.5_r8*(6._r8 * ncf)**(2._r8/3._r8)
         elseif( ncf.gt.1._r8/6._r8 .and. ncf.lt.1._r8 ) then
             phi = (acos(3._r8*(1._r8-ncf)/2._r8**(3._r8/2._r8))+4._r8*3.1415927_r8)/3._r8
             aist = (1._r8 - 4._r8 * cos(phi) * cos(phi))
         else
             aist = 1._r8
         endif
             aist = max(0._r8,min(aist,1._r8))
     elseif (iceopt.eq.5) then 
! set rh ice cloud fraction
             rhi= (qv+qi)/qsat * (esl/esi)
             rhdif= (rhi-rhmini) / (rhmaxi - rhmini)
             aist = min(1.0_r8, max(rhdif,0._r8)**2)

! limiter to remove empty cloud and ice with no cloud
! and set icecld fraction to mincld if ice exists

             if (qi.lt.minice) then
                aist=0._r8
             else
                aist=max(mincld,aist)
             endif

! enforce limits on icimr
             if (qi.ge.minice) then
                icimr=qi/aist

!minimum
                if (icimr.lt.qist_min) then
                   aist = max(0._r8,min(1._r8,qi/qist_min))
                endif
!maximum
                if (icimr.gt.qist_max) then
                   aist = max(0._r8,min(1._r8,qi/qist_max))
                endif

             endif
     endif 

   ! 0.999_r8 is added to prevent infinite 'ql_st' at the end of instratus_condensate
   ! computed after updating 'qi_st'.  

     aist = max(0._r8,min(aist,0.999_r8))

     aist_out(i) = aist

     enddo

   return
   end subroutine aist_vector

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

subroutine findsp_water_onelayer (lchnk, ncol, q, t, p, tsp, qsp)
!-----------------------------------------------------------------------
!
! Purpose:
!     find the wet bulb temperature for a given t and q
!     in a longitude height section
!     wet bulb temp is the temperature and spec humidity that is
!     just saturated and has the same enthalpy
!     if q > qs(t) then tsp > t and qsp = qs(tsp) < q
!     if q < qs(t) then tsp < t and qsp = qs(tsp) > q
!
! Method:
! a Newton method is used
! first guess uses an algorithm provided by John Petch from the UKMO
! we exclude points where the physical situation is unrealistic
! e.g. where the temperature is outside the range of validity for the
!      saturation vapor pressure, or where the water vapor pressure
!      exceeds the ambient pressure, or the saturation specific humidity is
!      unrealistic
!
! Author: P. Rasch
!
!-----------------------------------------------------------------------
!
!     input arguments
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns

   real(r8), intent(in) :: q(pcols)        ! water vapor (kg/kg)
   real(r8), intent(in) :: t(pcols)        ! temperature (K)
   real(r8), intent(in) :: p(pcols)        ! pressure    (Pa)
!
! output arguments
!
   real(r8), intent(out) :: tsp(pcols)      ! saturation temp (K)
   real(r8), intent(out) :: qsp(pcols)      ! saturation mixing ratio (kg/kg)
!
! local variables
!
   integer i                 ! work variable
!  integer k                 ! work variable
   logical lflg              ! work variable
   integer iter              ! work variable
   integer l                 ! work variable
   logical :: error_found

   real(r8) omeps                ! 1 minus epsilon
   real(r8) trinv                ! work variable
   real(r8) es                   ! sat. vapor pressure
   real(r8) desdt                ! change in sat vap pressure wrt temperature
!     real(r8) desdp                ! change in sat vap pressure wrt pressure
   real(r8) dqsdt                ! change in sat spec. hum. wrt temperature
   real(r8) dgdt                 ! work variable
   real(r8) g                    ! work variable
   real(r8) weight(pcols)        ! work variable
   real(r8) hlatsb               ! (sublimation)
   real(r8) hlatvp               ! (vaporization)
   real(r8) hltalt(pcols)   ! lat. heat. of vap.
   real(r8) tterm                ! work var.
   real(r8) qs                   ! spec. hum. of water vapor
   real(r8) tc                   ! crit temp of transition to ice
   real(r8) tt0                  ! Freezing temperature. Needed for findsp 

! work variables
   real(r8) t1, q1, dt, dq
   real(r8) dtm, dqm
   real(r8) qvd, a1, tmp
   real(r8) rair
   real(r8) r1b, c1, c2, c3
   real(r8) denom
   real(r8) dttol
   real(r8) dqtol
   integer doit(pcols)
   real(r8) enin(pcols), enout(pcols)
   real(r8) tlim(pcols)

   omeps = 1.0_r8 - epsqs
!  trinv = 1.0_r8/ttrice  ! put just before it is used, in case ttrice = 0.0  !++xl
   a1 = 7.5_r8*log(10._r8)
   rair =  287.04_r8
   c3 = rair*a1/cp
   dtm = 0._r8    ! needed for iter=0 blowup with f90 -ei
   dqm = 0._r8    ! needed for iter=0 blowup with f90 -ei
   dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
   dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
!  tmin = 173.16 ! the coldest temperature we can deal with
   tt0  = 273.15_r8
!
! max number of times to iterate the calculation
   iter = 8
!

!  Sungsu comment out k loop
!  do k = 1,pver

!
! first guess on the wet bulb temperature
!
      do i = 1,ncol

#ifdef DEBUG
         if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
            write (iulog,*) ' '
            write (iulog,*) ' level, t, q, p', t(i), q(i), p(i)
         endif
#endif
! limit the temperature range to that relevant to the sat vap pres tables
#if ( ! defined WACCM_MOZART )
         tlim(i) = min(max(t(i),173._r8),373._r8)
#else
         tlim(i) = min(max(t(i),128._r8),373._r8)
#endif
!        es = estblf(tlim(i))
         es = polysvp(tlim(i),0)  !++xl
         denom = p(i) - omeps*es
         qs = epsqs*es/denom
         doit(i) = 0
         enout(i) = 1._r8
! make sure a meaningful calculation is possible
         if (p(i) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then
!
! Saturation specific humidity
!
             qs = min(epsqs*es/denom,1._r8)
!
! "generalized" analytic expression for t derivative of es
!  accurate to within 1 percent for 173.16 < t < 373.16

!++xl
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
!            tc     = tlim(i) - tt0
!            lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
!            weight(i) = min(-tc*trinv,1.0_r8)
!            hlatsb = hlatv + weight(i)*hlatf
!            hlatvp = hlatv - 2369.0_r8*tc
!            if (tlim(i) < tt0) then
!               hltalt(i) = hlatsb
!            else
!               hltalt(i) = hlatvp
!            end if
!
! No icephs or water to ice transition
!
             hlatvp = hlatv - 2369.0*(tlim(i)-tt0)
             hlatsb = hlatv
             if (tlim(i) < tt0) then
               hltalt(i) = hlatsb
             else
               hltalt(i) = hlatvp
             end if
!--xl
             enin(i) = cp*tlim(i) + hltalt(i)*q(i)

! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
             tmp =  q(i) - qs
             c1 = hltalt(i)*c3
             c2 = (tlim(i) + 36._r8)**2
             r1b    = c2/(c2 + c1*qs)
             qvd   = r1b*tmp
             tsp(i) = tlim(i) + ((hltalt(i)/cp)*qvd)
#ifdef DEBUG
             if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
                write (iulog,*) ' relative humidity ', q(i)/qs
                write (iulog,*) ' first guess ', tsp(i)
             endif
#endif
!            es = estblf(tsp(i))
             es = polysvp(tsp(i),0)  !++xl
             qsp(i) = min(epsqs*es/(p(i) - omeps*es),1._r8)
          else
             doit(i) = 1
             tsp(i) = tlim(i)
             qsp(i) = q(i)
             enin(i) = 1._r8
          endif
       end do   ! end do i
!
! now iterate on first guess
!
      do l = 1, iter
         dtm = 0
         dqm = 0
         do i = 1,ncol
            if (doit(i) == 0) then
!              es = estblf(tsp(i))
               es = polysvp(tsp(i),0)  !++xl
!
! Saturation specific humidity
!
               qs = min(epsqs*es/(p(i) - omeps*es),1._r8)
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!
!++xl
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
!              tc     = tsp(i) - tt0
!              lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
!              weight(i) = min(-tc*trinv,1.0_r8)
!              hlatsb = hlatv + weight(i)*hlatf
!              hlatvp = hlatv - 2369.0_r8*tc
!              if (tsp(i) < tt0) then
!                 hltalt(i) = hlatsb
!              else
!                 hltalt(i) = hlatvp
!              end if
!              if (lflg) then
!                 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5))))
!              else
!                 tterm = 0.0_r8
!              end if
!              desdt = hltalt(i)*es/(rgasv*tsp(i)*tsp(i)) + tterm*trinv
!
! No icephs or water to ice transition
!
               hlatvp = hlatv - 2369.0*(tsp(i)-tt0)
               hlatsb = hlatv
               if (tsp(i) < tt0) then
                 hltalt(i) = hlatsb
               else
                 hltalt(i) = hlatvp
               end if
               desdt = hltalt(i)*es/(rgasv*tsp(i)*tsp(i))
!--xl
               dqsdt = (epsqs + omeps*qs)/(p(i) - omeps*es)*desdt
!              g = cp*(tlim(i)-tsp(i)) + hltalt(i)*q(i)- hltalt(i)*qsp(i)
               g = enin(i) - (cp*tsp(i) + hltalt(i)*qsp(i))
               dgdt = -(cp + hltalt(i)*dqsdt)
               t1 = tsp(i) - g/dgdt
               dt = abs(t1 - tsp(i))/t1
               tsp(i) = max(t1,tmin)
!              es = estblf(tsp(i))
               es = polysvp(tsp(i),0)  !++xl
               q1 = min(epsqs*es/(p(i) - omeps*es),1._r8)
               dq = abs(q1 - qsp(i))/max(q1,1.e-12_r8)
               qsp(i) = q1
#ifdef DEBUG
               if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
                  write (iulog,*) ' rel chg lev, iter, t, q ', l, dt, dq, g
               endif
#endif
               dtm = max(dtm,dt)
               dqm = max(dqm,dq)
! if converged at this point, exclude it from more iterations
               if (dt < dttol .and. dq < dqtol) then
                  doit(i) = 2
               endif
               enout(i) = cp*tsp(i) + hltalt(i)*qsp(i)
! bail out if we are too near the end of temp range
#if ( ! defined WACCM_MOZART )
               if (tsp(i) < 174.16_r8) then
#else
               if (tsp(i) < 130.16_r8) then
#endif
                  doit(i) = 4
               endif
            else
            endif
         end do              ! do i = 1,ncol

         if (dtm < dttol .and. dqm < dqtol) then
            go to 10
         endif

      end do                 ! do l = 1,iter
10    continue

      error_found = .false.
      if (dtm > dttol .or. dqm > dqtol) then
         do i = 1,ncol
            if (doit(i) == 0) error_found = .true.
         end do
         if (error_found) then
            do i = 1,ncol
               if (doit(i) == 0) then
                  write (iulog,*) ' findsp not converging at point i ', i
                  write (iulog,*) ' t, q, p, enin ', t(i), q(i), p(i), enin(i)
                  write (iulog,*) ' tsp, qsp, enout ', tsp(i), qsp(i), enout(i)
                  call endrun ('FINDSP')
               endif
            end do
         endif
      endif
      do i = 1,ncol
         if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
            error_found = .true.
         endif
      end do
      if (error_found) then
         do i = 1,ncol
            if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
               write (iulog,*) ' the enthalpy is not conserved for point ', &
                  i, enin(i), enout(i)
               write (iulog,*) ' t, q, p, enin ', t(i), q(i), p(i), enin(i)
               write (iulog,*) ' tsp, qsp, enout ', tsp(i), qsp(i), enout(i)
               call endrun ('FINDSP')
            endif
         end do
      endif
      
 ! end do                    ! level loop (k=1,pver)

   return
end subroutine findsp_water_onelayer
!--xl

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

subroutine findsp_water (lchnk, ncol, q, t, p, tsp, qsp)
!-----------------------------------------------------------------------
!
! Purpose:
!     find the wet bulb temperature for a given t and q
!     in a longitude height section
!     wet bulb temp is the temperature and spec humidity that is
!     just saturated and has the same enthalpy
!     if q > qs(t) then tsp > t and qsp = qs(tsp) < q
!     if q < qs(t) then tsp < t and qsp = qs(tsp) > q
!
! Method:
! a Newton method is used
! first guess uses an algorithm provided by John Petch from the UKMO
! we exclude points where the physical situation is unrealistic
! e.g. where the temperature is outside the range of validity for the
!      saturation vapor pressure, or where the water vapor pressure
!      exceeds the ambient pressure, or the saturation specific humidity is
!      unrealistic
!
! Author: P. Rasch
!
!-----------------------------------------------------------------------
!
!     input arguments
!
   integer, intent(in) :: lchnk                 ! chunk identifier
   integer, intent(in) :: ncol                  ! number of atmospheric columns

   real(r8), intent(in) :: q(pcols,pver)        ! water vapor (kg/kg)
   real(r8), intent(in) :: t(pcols,pver)        ! temperature (K)
   real(r8), intent(in) :: p(pcols,pver)        ! pressure    (Pa)
!
! output arguments
!
   real(r8), intent(out) :: tsp(pcols,pver)      ! saturation temp (K)
   real(r8), intent(out) :: qsp(pcols,pver)      ! saturation mixing ratio (kg/kg)
!
! local variables
!
   integer i                 ! work variable
   integer k                 ! work variable
   logical lflg              ! work variable
   integer iter              ! work variable
   integer l                 ! work variable
   logical :: error_found

   real(r8) omeps                ! 1 minus epsilon
   real(r8) trinv                ! work variable
   real(r8) es                   ! sat. vapor pressure
   real(r8) desdt                ! change in sat vap pressure wrt temperature
!     real(r8) desdp                ! change in sat vap pressure wrt pressure
   real(r8) dqsdt                ! change in sat spec. hum. wrt temperature
   real(r8) dgdt                 ! work variable
   real(r8) g                    ! work variable
   real(r8) weight(pcols)        ! work variable
   real(r8) hlatsb               ! (sublimation)
   real(r8) hlatvp               ! (vaporization)
   real(r8) hltalt(pcols,pver)   ! lat. heat. of vap.
   real(r8) tterm                ! work var.
   real(r8) qs                   ! spec. hum. of water vapor
   real(r8) tc                   ! crit temp of transition to ice
   real(r8) tt0                  ! Freezing temperature. Needed for findsp 

! work variables
   real(r8) t1, q1, dt, dq
   real(r8) dtm, dqm
   real(r8) qvd, a1, tmp
   real(r8) rair
   real(r8) r1b, c1, c2, c3
   real(r8) denom
   real(r8) dttol
   real(r8) dqtol
   integer doit(pcols)
   real(r8) enin(pcols), enout(pcols)
   real(r8) tlim(pcols)

   omeps = 1.0_r8 - epsqs
!  trinv = 1.0_r8/ttrice  ! put just before it is used, in case ttrice = 0.0  !++xl
   a1 = 7.5_r8*log(10._r8)
   rair =  287.04_r8
   c3 = rair*a1/cp
   dtm = 0._r8    ! needed for iter=0 blowup with f90 -ei
   dqm = 0._r8    ! needed for iter=0 blowup with f90 -ei
   dttol = 1.e-4_r8 ! the relative temp error tolerance required to quit the iteration
   dqtol = 1.e-4_r8 ! the relative moisture error tolerance required to quit the iteration
!  tmin = 173.16 ! the coldest temperature we can deal with
   tt0  = 273.15_r8 
!
! max number of times to iterate the calculation
   iter = 8
!

   do k = 1,pver

!
! first guess on the wet bulb temperature
!
      do i = 1,ncol

#ifdef DEBUG
         if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
            write (iulog,*) ' '
            write (iulog,*) ' level, t, q, p', t(i,k), q(i,k), p(i,k)
         endif
#endif
! limit the temperature range to that relevant to the sat vap pres tables
#if ( ! defined WACCM_MOZART )
         tlim(i) = min(max(t(i,k),173._r8),373._r8)
#else
         tlim(i) = min(max(t(i,k),128._r8),373._r8)
#endif
!        es = estblf(tlim(i))
         es = polysvp(tlim(i),0)  !++xl
         denom = p(i,k) - omeps*es
         qs = epsqs*es/denom
         doit(i) = 0
         enout(i) = 1._r8
! make sure a meaningful calculation is possible
         if (p(i,k) > 5._r8*es .and. qs > 0._r8 .and. qs < 0.5_r8) then
!
! Saturation specific humidity
!
             qs = min(epsqs*es/denom,1._r8)
!
! "generalized" analytic expression for t derivative of es
!  accurate to within 1 percent for 173.16 < t < 373.16

!++xl
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
!            tc     = tlim(i) - tt0
!            lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
!            weight(i) = min(-tc*trinv,1.0_r8)
!            hlatsb = hlatv + weight(i)*hlatf
!            hlatvp = hlatv - 2369.0_r8*tc
!            if (tlim(i) < tt0) then
!               hltalt(i,k) = hlatsb
!            else
!               hltalt(i,k) = hlatvp
!            end if
!
! No icephs or water to ice transition
!
             hlatvp = hlatv - 2369.0*(tlim(i)-tt0)
             hlatsb = hlatv
             if (tlim(i) < tt0) then
               hltalt(i,k) = hlatsb
             else
               hltalt(i,k) = hlatvp
             end if
!--xl
             enin(i) = cp*tlim(i) + hltalt(i,k)*q(i,k)

! make a guess at the wet bulb temp using a UKMO algorithm (from J. Petch)
             tmp =  q(i,k) - qs
             c1 = hltalt(i,k)*c3
             c2 = (tlim(i) + 36._r8)**2
             r1b    = c2/(c2 + c1*qs)
             qvd   = r1b*tmp
             tsp(i,k) = tlim(i) + ((hltalt(i,k)/cp)*qvd)
#ifdef DEBUG
             if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
                write (iulog,*) ' relative humidity ', q(i,k)/qs
                write (iulog,*) ' first guess ', tsp(i,k)
             endif
#endif
!            es = estblf(tsp(i,k))
             es = polysvp(tsp(i,k),0)  !++xl
             qsp(i,k) = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
          else
             doit(i) = 1
             tsp(i,k) = tlim(i)
             qsp(i,k) = q(i,k)
             enin(i) = 1._r8
          endif
       end do   ! end do i
!
! now iterate on first guess
!
      do l = 1, iter
         dtm = 0
         dqm = 0
         do i = 1,ncol
            if (doit(i) == 0) then
!              es = estblf(tsp(i,k))
               es = polysvp(tsp(i,k),0)  !++xl
!
! Saturation specific humidity
!
               qs = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
!
! "generalized" analytic expression for t derivative of es
! accurate to within 1 percent for 173.16 < t < 373.16
!
!++xl
! Weighting of hlat accounts for transition from water to ice
! polynomial expression approximates difference between es over
! water and es over ice from 0 to -ttrice (C) (min of ttrice is
! -40): required for accurate estimate of es derivative in transition
! range from ice to water also accounting for change of hlatv with t
! above freezing where const slope is given by -2369 j/(kg c) = cpv - cw
!
!              tc     = tsp(i,k) - tt0
!              lflg   = (tc >= -ttrice .and. tc < 0.0_r8)
!              weight(i) = min(-tc*trinv,1.0_r8)
!              hlatsb = hlatv + weight(i)*hlatf
!              hlatvp = hlatv - 2369.0_r8*tc
!              if (tsp(i,k) < tt0) then
!                 hltalt(i,k) = hlatsb
!              else
!                 hltalt(i,k) = hlatvp
!              end if
!              if (lflg) then
!                 tterm = pcf(1) + tc*(pcf(2) + tc*(pcf(3)+tc*(pcf(4) + tc*pcf(5))))
!              else
!                 tterm = 0.0_r8
!              end if
!              desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k)) + tterm*trinv
!
! No icephs or water to ice transition
!
               hlatvp = hlatv - 2369.0*(tsp(i,k)-tt0)
               hlatsb = hlatv
               if (tsp(i,k) < tt0) then
                 hltalt(i,k) = hlatsb
               else
                 hltalt(i,k) = hlatvp
               end if
               desdt = hltalt(i,k)*es/(rgasv*tsp(i,k)*tsp(i,k))
!--xl
               dqsdt = (epsqs + omeps*qs)/(p(i,k) - omeps*es)*desdt
!              g = cp*(tlim(i)-tsp(i,k)) + hltalt(i,k)*q(i,k)- hltalt(i,k)*qsp(i,k)
               g = enin(i) - (cp*tsp(i,k) + hltalt(i,k)*qsp(i,k))
               dgdt = -(cp + hltalt(i,k)*dqsdt)
               t1 = tsp(i,k) - g/dgdt
               dt = abs(t1 - tsp(i,k))/t1
               tsp(i,k) = max(t1,tmin)
!              es = estblf(tsp(i,k))
               es = polysvp(tsp(i,k),0)  !++xl
               q1 = min(epsqs*es/(p(i,k) - omeps*es),1._r8)
               dq = abs(q1 - qsp(i,k))/max(q1,1.e-12_r8)
               qsp(i,k) = q1
#ifdef DEBUG
               if ( (lchnk == lchnklook(nlook) ) .and. (i == icollook(nlook) ) ) then
                  write (iulog,*) ' rel chg lev, iter, t, q ', l, dt, dq, g
               endif
#endif
               dtm = max(dtm,dt)
               dqm = max(dqm,dq)
! if converged at this point, exclude it from more iterations
               if (dt < dttol .and. dq < dqtol) then
                  doit(i) = 2
               endif
               enout(i) = cp*tsp(i,k) + hltalt(i,k)*qsp(i,k)
! bail out if we are too near the end of temp range
#if ( ! defined WACCM_MOZART )
               if (tsp(i,k) < 174.16_r8) then
#else
               if (tsp(i,k) < 130.16_r8) then
#endif
                  doit(i) = 4
               endif
            else
            endif
         end do              ! do i = 1,ncol

         if (dtm < dttol .and. dqm < dqtol) then
            go to 10
         endif

      end do                 ! do l = 1,iter
10    continue

      error_found = .false.
      if (dtm > dttol .or. dqm > dqtol) then
         do i = 1,ncol
            if (doit(i) == 0) error_found = .true.
         end do
         if (error_found) then
            do i = 1,ncol
               if (doit(i) == 0) then
                  write (iulog,*) ' findsp not converging at point i ', i
                  write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
                  write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
                  call endrun ('FINDSP')
               endif
            end do
         endif
      endif
      do i = 1,ncol
         if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
            error_found = .true.
         endif
      end do
      if (error_found) then
         do i = 1,ncol
            if (doit(i) == 2 .and. abs((enin(i)-enout(i))/(enin(i)+enout(i))) > 1.e-4_r8) then
               write (iulog,*) ' the enthalpy is not conserved for point ', &
                  i, enin(i), enout(i)
               write (iulog,*) ' t, q, p, enin ', t(i,k), q(i,k), p(i,k), enin(i)
               write (iulog,*) ' tsp, qsp, enout ', tsp(i,k), qsp(i,k), enout(i)
               call endrun ('FINDSP')
            endif
         end do
      endif
      
   end do                    ! level loop (k=1,pver)

   return
end subroutine findsp_water
!--xl

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

      SUBROUTINE gaussj(a,n,np,b,m,mp)
      INTEGER m,mp,n,np,NMAX
      real(r8) a(np,np),b(np,mp)
      real(r8) aa(np,np),bb(np,mp)
      PARAMETER (NMAX=50)
      INTEGER i,icol,irow,j,k,l,ll,ii,jj,indxc(NMAX),indxr(NMAX),ipiv(NMAX)
      real(r8) big,dum,pivinv

      aa(:,:) = a(:,:)
      bb(:,:) = b(:,:)

      do 11 j=1,n
        ipiv(j)=0
11    continue
      do 22 i=1,n
        big=0.
        do 13 j=1,n
          if(ipiv(j).ne.1)then
            do 12 k=1,n
              if (ipiv(k).eq.0) then
                if (abs(a(j,k)).ge.big)then
                  big=abs(a(j,k))
                  irow=j
                  icol=k
                endif
              else if (ipiv(k).gt.1) then
                write(iulog,*) 'singular matrix in gaussj 1'
                do ii = 1, np
                do jj = 1, np
                   write(iulog,*) ii, jj, aa(ii,jj), bb(ii,1)
                end do
                end do   
                call endrun
              endif
12          continue
          endif
13      continue
        ipiv(icol)=ipiv(icol)+1
        if (irow.ne.icol) then
          do 14 l=1,n
            dum=a(irow,l)
            a(irow,l)=a(icol,l)
            a(icol,l)=dum
14        continue
          do 15 l=1,m
            dum=b(irow,l)
            b(irow,l)=b(icol,l)
            b(icol,l)=dum
15        continue
        endif
        indxr(i)=irow
        indxc(i)=icol
        if (a(icol,icol).eq.0.) then
            write(iulog,*) 'singular matrix in gaussj 2'
            do ii = 1, np
            do jj = 1, np
               write(iulog,*) ii, jj, aa(ii,jj), bb(ii,1)
            end do
            end do   
            call endrun
        endif 
        pivinv=1./a(icol,icol)
        a(icol,icol)=1.
        do 16 l=1,n
          a(icol,l)=a(icol,l)*pivinv
16      continue
        do 17 l=1,m
          b(icol,l)=b(icol,l)*pivinv
17      continue
        do 21 ll=1,n
          if(ll.ne.icol)then
            dum=a(ll,icol)
            a(ll,icol)=0.
            do 18 l=1,n
              a(ll,l)=a(ll,l)-a(icol,l)*dum
18          continue
            do 19 l=1,m
              b(ll,l)=b(ll,l)-b(icol,l)*dum
19          continue
          endif
21      continue
22    continue
      do 24 l=n,1,-1
        if(indxr(l).ne.indxc(l))then
          do 23 k=1,n
            dum=a(k,indxr(l))
            a(k,indxr(l))=a(k,indxc(l))
            a(k,indxc(l))=dum
23        continue
        endif
24    continue

      return
      end subroutine gaussj

   ! ----------------- !
   ! End of subroutine !
   ! ----------------- !

end module cldwat2m_macro
